Enzyme utilization and saturation for R. eutropha

1 Description

This R notebook is a bioinformatics pipeline to analyze protein saturation/under-utilization with a resource allocation model for the chemolithoautotroph Ralstonia eutropha (a.k.a. Cupriavidus necator).

2 Libraries

suppressPackageStartupMessages({
  library(lattice)
  library(latticeExtra)
  library(latticetools)
  library(tidyverse)
  library(stringi)
})

3 Data import

Define the data source directories. Some of them are external in the sense of not included in the accompanying data folder of this R notebook. These are located in the accompanying github repository for the resource allocation model that was used here. The resource allocation model can be found at my fork of Bacterial-RBA-models.

Reutropha_proteomics <- "../data/input/Ralstonia_eutropha.Rdata"
model_reactions <- "../data/input/model_reactions.csv"
simulation_dir <- "../data/simulation/substrate_limitation/"
source("read_rba_result.R")

Read simulation data.

# read simulation results
df_flux <- read_rba_result(list.files(simulation_dir, pattern = "fluxes_.*.tsv$", full.names = TRUE))
df_prot <- read_rba_result(list.files(simulation_dir, pattern = "proteins_.*.tsv", full.names = TRUE))
df_macr <- read_rba_result(list.files(simulation_dir, pattern = "macroprocesses_.*.tsv", full.names = TRUE))

4 Overview on substrate uptake, growth, and yield

After running a set of simulations in RBApy that simulate increasing substrate limitation, we can plot the substrate uptake rate q, yield Y in gram biomass per gram substrate, and growth rate µ. Unlike genome scale models, growth becomes limited by the maximum amount of proteins that a cell can synthesize. If cells would not be protein-limited, or proteins would catalyze reactions infinitely fast, no such limitation would take place and growth rate would scale linearly with substrate concentration. This is the situation in FBA simulation.

# rearrange some rows (mu, qS) to columns
df_macr <- df_macr %>% filter(!grepl("test_process", key)) %>% 
  spread(key, value) %>%
  
  # add type of simulation
  mutate(substrate = case_when(
    carbon_source == "for" ~ "formate",
    carbon_source == "succ" ~ "succinate",
    carbon_source == "fru" & nitrogen_conc == 18.7 ~ "fructose",
    carbon_source == "fru" & nitrogen_conc != 18.7 ~ "ammonium"
  )) %>%
  
  # add uptake rate in g/gDCW*h instead of mmol
  mutate(qS_g_gDCW_h = case_when(
    substrate == "formate" ~ qS*0.04603,
    substrate == "succinate" ~ qS*0.11809,
    substrate == "fructose" ~ qS*0.18016,
    substrate == "ammonium" ~ qS*0.05349
  ))

First we can have a look at how growth rate levels off with increasing substrate concentration in mmol/L. Note: this is not equal to substrate uptake rate. Substrate uptake rate and growth rate should have an almost linear relationship. We can log-transform the carbon concentration and fit a linear model to predict the substrate concentration required to obtain a certain substrate uptake rate and growth rate.

# copy nitrogen to carbon concentration for ammonium limitation,
# just for plotting purposes
df_macr_viz <- df_macr %>%
  mutate(carbon_conc = case_when(
    substrate == "ammonium" ~ nitrogen_conc,
    TRUE ~ carbon_conc
  ))

plot_mu_qs_lin <- xyplot(mu ~ carbon_conc | substrate, df_macr_viz,
    par.settings = custom.colorblind(),
    between = list(x = 0.5, y = 0.5),
    layout = c(4,1), lwd = 1.5, pch = 19,
    xlab = expression("S [mM]"),
    ylab = expression('µ [h'^'-1'*']'),
    scales = list(alternating = FALSE),
    panel = function(x, y, ...) {
      panel.grid(h = -1, v = -1, col = grey(0.9))
      panel.xyplot(x, y, cex = 0.9, ...)
      #panel.lmlineq(x[1:4], y[1:4], fontfamily = "FreeSans", ...)
    }
  )

plot_mu_qs_log <- xyplot(log10(carbon_conc) ~ qS | substrate,
    df_macr_viz,
    par.settings = custom.colorblind(),
    between = list(x = 0.5, y = 0.5), #xlim = c(0, 5),
    layout = c(4,1), lwd = 1.5, pch = 19,
    xlab = expression('q'[S]*' mmol g DCW'^-1*'h'^-1),
    ylab =  expression('log'[10]*' S [mM]'),
    scales = list(alternating = FALSE),
    panel = function(x, y, ...) {
      panel.grid(h = -1, v = -1, col = grey(0.9))
      panel.xyplot(x, y, cex = 0.9, ...)
      panel.lmlineq(x[1:6], y[1:6], fontfamily = "FreeSans", ...)
    }
  )

print(plot_mu_qs_lin, split = c(1,1,1,2), more = TRUE)
print(plot_mu_qs_log, split = c(1,2,1,2))

df_macr_viz %>%
  group_by(substrate) %>%
  # fit linear model to substrate uptake rate vs concentration
  summarize(
    slope = lm(x ~ y, data = list(x = log10(carbon_conc), y = qS))$coeff[2],
    offset = lm(x ~ y, data = list(x = log10(carbon_conc), y = qS))$coeff[1]
  ) %>% mutate(model = paste0("c = 10^(", round(offset, 3), " + ", round(slope, 3), "*qS)"))

Create a Herbert-Pirt plot for each condition (growth rate versus substrate uptake rate). This plot would show a change in yield by a ‘kink’ of the data points.

xyplot(qS_g_gDCW_h ~ mu | substrate, df_macr,
  par.settings = custom.colorblind(),
  between = list(x = 0.5, y = 0.5),
  layout = c(4,1), lwd = 1.5, pch = 19,
  ylab = expression("q"[S]*" [g h"^-1*" gDCW"^-1*"]"),
  xlab = expression('µ [h'^'-1'*']'),
  scales = list(alternating = FALSE),
  panel = function(x, y, ...) {
    panel.grid(h = -1, v = -1, col = grey(0.9))
    panel.xyplot(x, y, cex = 0.9, ...)
    # displaying maintenance and yield coefficients
    coef <- lm(y ~ x, data.frame(x, y))$coeff
    panel.text(median(x), 2.7, 
      paste("ms =", round(coef[[1]], 3), "g h-1 g_DCW-1"), 
      col = grey(0.3), cex = 0.7)
    panel.text(median(x), 2.4, paste(expression("Yx/S ="), 
        round(1/coef[[2]], 3), "g_DCW g_S-1"), 
      col = grey(0.3), cex = 0.7)
  }
)

5 Resource allocation in terms of protein mass

To determine the true allocation of protein resources per compartment, but also the true cost of protein per process, we need to convert the predicted concentration of proteins in mmol per gDCW to g per gDCW, simply by multiplying protein concentration with the molecular weight of a protein (g/mol, converted to g/mmol). We can then also easily transform g/gDCW to mass fraction by dividing individual protein concentrations by the sum of all protein concentrations. The protein mass fraction is dimensionless. The only parameter required for this transformation is the molecular weight per protein which is available from uniprot. We can for example take the protein annotation table that is automatically downloaded during RBApy model generation.

# import downloaded Ralstonia protein annotation from uniprot
df_uniprot <- read_tsv("../data/input/uniprot.csv", col_types = cols()) %>%
  mutate(locus_tag = stri_extract_first(`Gene names`, regex = "H16_[AB][0-9]{4}|PHG[0-9]{3}"))

# merge predicted protein allocation with molecular weight info from uniprot
df_prot <- left_join(df_prot, select(df_uniprot, locus_tag, Length, Mass),
  by = c("key" = "locus_tag")) %>%
  
  # calculate predicted protein mass in g/gDCW using MW in g/mmol, and mass fraction
  group_by(simulation) %>% mutate(
    predicted_mass_g_gDCW = value * Mass / 1000)

# test if mass fractions sum to reasonable value
df_prot %>% summarize(
  predicted_mass_g_gDCW = sum(predicted_mass_g_gDCW, na.rm = TRUE))

---

The simulated protein mass per gDCW changes with growth rate, and it is considerably lower then the estimated ~0.65-0.68 g total protein/gDCW that was previously estimated for bacteria (see e.g. Park et al., biomass composition for Ralstonia eutropha model, or Touloupakis et al., biomass composition for cyanobacteria). One reason is that above numbers only include enzymatic proteins (the ones represented by the GSM). For machinery such as ribosomal proteins, a separate calculation needs to be performed. The model returns estimated concentration of machineries for replication, transcription, translation, and protein folding. The associated proteins can be imported from the model folder and the total mass estimated using molecular weight and subunit stoichiometry as done before for enzymatic proteins.

machinery_names <- c("replication", "transcription", "ribosome", "chaperones")
machinery_tables <- paste0("../data/simulation/macro_machines/", machinery_names, ".tsv")

df_macr <- df_macr %>% 
  
  # remove unused columns and rename important ones
  select(-P_ENZ, -P_RNADEG) %>%
  rename(replication = P_REP, transcription = P_TSC, 
    ribosome = P_TA, chaperones = P_CHP, substrate_uptake_rate = qS,
    predicted_growth_rate = mu
  ) %>%
  
  # gather individual machineries in one column
  gather(key = machine, value = predicted_mass_mmol_gDCW, 
    chaperones:transcription) %>%
  
  # round uptake rates for merging
  mutate(substrate_uptake_rate = round(substrate_uptake_rate, 3))
  

df_machinery <- lapply(machinery_tables, read_tsv, col_type = cols()) %>% 
  bind_rows(.id = "machine") %>%
  mutate(machine = recode(machine, !!!setNames(machinery_names, 1:4))) %>%
  select(-`Entry name`, -Sequence, -Cofactor, -`EC number`, -`Organism ID`, `Organism`,
    -`Catalytic activity`, -Status) %>%
  
  # join with prediction of molecular machine concentration (mmol/gDCW)
  left_join(df_macr, by = "machine") %>%
  
  # calculate predicted protein mass in g/gDCW using MW in g/mmol, and mass fraction
  group_by(simulation) %>% mutate(
    predicted_mass_g_gDCW = predicted_mass_mmol_gDCW * Stoichiometry * Mass / 1000
  )

# test if mass fractions sum to reasonable value
df_machinery %>% 
  summarize(sum_g_gDCW = sum(predicted_mass_g_gDCW, na.rm = TRUE)) %>%
  ungroup %>% head(10)

If the utilized protein for enzymes and machinery is summed up, it does not exceed ~0.25 g/gDCW. The reason for this is that up to 60% of the protein mass is not modeled (non enzymatic, NE proteome, 57.7% at µ = 0, see R notebook Ralstonia-model-constraints), and only around 68% of the total cell mass is protein. If we consider this, than total protein mass can be estimated as m = 0.25/(1-0.6) = 0.625 g/gDCW. This is much closer to the estimated 0.68 g protein/gDCW. The difference can be attributed to default values of amino acid concentration in RBApy that determine the size of the ‘protein pool’. It is important to note that the sum of utilizable proteins is not a constant value but a linear function of growth rate. The total proteome pool including non-enzymatic proteins is, however, constant.

6 Correlation between predicted and experimentally determined proteome

To compare the predicted and experimental proteome composition, we load the required proteomics data. Proteomics data are mass spectrometry measurements with label-free quantification of peptides. Protein quantification was performed by summing up all peptide intensities per annotated protein. The proteomic measurement unit, mass fraction, can be easily transformed to g/gDCW by multiplying mass fraction with the total protein mass (0.68 g/gDCW) used in RBApy simulations. Or vice versa by converting RBApy protein concentration (mmol/gDCW) to mass fraction.

To allow a fair comparison between measured and predicted data, it is necessary to aggregate (e.g. sum up) all protein abundances allocated to one reaction. The reason is that the model will only predict protein abundance of the first of a range of iso-enzymes for a particular reaction, while in reality another iso-enzyme might be more abundant (carry the majority of flux). This would lead to lower correlation between measured and predicted protein concentrations. This is not necessary for machinery proteins, they have no iso-enzymes and are used under all conditions.

6.1 Combine simulations and experimental data

The proteomics data must be merged with RBApy simulation results using matching conditions. First, proteomics data are loaded and prepared for merging.

Step 1: load proteomics data

load(Reutropha_proteomics)

# pick a condition matching simulations
Ralstonia_eutropha <- Ralstonia_eutropha %>%
  
  # round substrate uptake rate
  mutate(substrate_uptake_rate = round(substrate_uptake_rate, 3)) %>%
  
  # select only required columns
  rename(growth_rate = growthrate) %>%
  select(uniprot, locus_tag, protein, condition, substrate, substrate_uptake_rate,
    growth_rate, COG_Process, R1:R4) %>%
  
  # turn raw intensity measurements into mass in g per gDCW (assuming a 
  # total protein concentration of 0.68 g/gDCW)
  group_by(condition) %>%
  mutate(across(matches("R[1234]"), function(x) x/sum(x, na.rm = TRUE)*0.68)) %>%
  gather(key = "replicate", value = "mass_g_gDCW", R1:R4)

Step 2: Load gene reaction associations obtained from genome scale model

df_model_reactions <- read_csv(model_reactions, col_types = cols()) %>%
  
  # filter for reactions with gene associations
  select(reaction_id, reaction_name, genes, groups) %>% separate_rows(genes, sep = ", ") %>%
  filter(!is.na(genes)) %>% 
  rename(model_group = groups) %>%
  
  # add a more general groups description
  mutate(model_group_basic = case_when(
    grepl("Phenylala|Valine|Tyrosine|Glutamate|Glycine|Tryptophan|Methionine|
          Cysteine|Alanine|Histidine|Arginine|Lysine", model_group) ~ "Amino acid",
    grepl("Pentose|Calvin", model_group) ~ "PPP + Calvin cycle",
    model_group == "Citric Acid Cycle" ~ "Citric Acid Cycle",
    model_group == "Glycolysis/Gluconeogenesis" ~ "Glycolysis/Gluconeogenesis",
    model_group == "Glyoxylate and Dicarboxylate metabolism" ~ "Autotrophic energy",
    model_group == "Oxidative Phosphorylation" ~ "Oxidative Phosphorylation",
    TRUE ~ "Other"
  ))

Step 3: Select and rename conditions from RBA simulation

# add type of substrate limitation
add_cond <- function(df) {
  df %>% mutate(substrate = case_when(
    carbon_source == "succ" ~ "succinate",
    carbon_source == "for" ~ "formate",
    carbon_source == "fru" & nitrogen_conc == 18.7 ~ "fructose",
    TRUE ~ "ammonium",
  ))
}

# add substrate uptake rate
df_substrate_uptake <- df_macr %>% 
  select(simulation, substrate_uptake_rate) %>% distinct
df_prot <- df_prot %>% add_cond %>% left_join(df_substrate_uptake, by = "simulation")
df_flux <- df_flux %>% add_cond %>% left_join(df_substrate_uptake, by = "simulation")

Step 4: Merge protein measurements and predictions into master tables

The first step is to merge the tables for machinery proteins, that means proteins related to replication, transcription, translation, and protein folding. These don’t require allocation of protein mass to reactions, and merging becomes simply an operation on enzyme IDs and conditions.

# join with proteomics data
df_machinery <- df_machinery %>%
  left_join(Ralstonia_eutropha,
    by = c("Entry" = "uniprot", "substrate", "substrate_uptake_rate"))

The second table for all enzymatic proteins requires the allocation of estimated protein mass to enzymes. One option for the future is to retrieve these values directly from RBApy, but this is not implemented yet.

df_prot_comp <- df_model_reactions %>%
  
  # join with proteomics data
  left_join(Ralstonia_eutropha, by = c("genes" = "locus_tag")) %>%
  complete(nesting(genes, reaction_id, reaction_name, model_group, model_group_basic),
    nesting(condition, substrate, substrate_uptake_rate, growth_rate, replicate)) %>%
  filter(!is.na(substrate)) %>%
  
  # join with simulation data
  left_join(df_prot, by = c("genes" = "key", "substrate", "substrate_uptake_rate")) %>%
  
  # determine number of reactions per protein
  group_by(condition, genes, replicate) %>% 
  mutate(n_reactions = length(reaction_id)) %>%
  
  # calculate protein mass in g/gDCW
  ungroup %>% mutate(
    predicted_mass_g_gDCW = predicted_mass_g_gDCW/n_reactions,
    mass_g_gDCW = mass_g_gDCW/n_reactions
  ) %>%
  
  # summarize by summing up protein abundance per reaction (NA treated as zero)
  group_by(condition, reaction_id, reaction_name, model_group,
    model_group_basic, substrate, substrate_uptake_rate, growth_rate, replicate) %>% 
  summarize(
    predicted_mass_g_gDCW = sum(predicted_mass_g_gDCW, na.rm = TRUE),
    mass_g_gDCW = sum(mass_g_gDCW, na.rm = TRUE)
  ) %>%
  filter(mass_g_gDCW != 0) %>%
  
  # add predicted growth rate to experimental
  left_join(
    df_machinery %>% ungroup %>%
    select(substrate, substrate_uptake_rate, predicted_growth_rate) %>%
    filter(!duplicated(substrate_uptake_rate))
  ) %>%
  
  # add predicted fluxes per reaction and condition
  left_join(
    df_flux %>% ungroup %>%
    select(key, value, substrate, substrate_uptake_rate) %>%
    rename(reaction_id = key, flux_mmol_gDCW_h = value)
  )

`summarise()` has grouped output by 'condition', 'reaction_id', 'reaction_name', 'model_group', 'model_group_basic', 'substrate', 'substrate_uptake_rate', 'growth_rate'. You can override using the `.groups` argument.
Joining, by = c("substrate", "substrate_uptake_rate")
Joining, by = c("reaction_id", "substrate", "substrate_uptake_rate")

Now we perform a test. We check if all protein abundances allocated to reactions sum to a reasonable value as we would expect. This value would be the total enzymatic protein mass in g/gDCW, per condition and replicate, and could reach up to 0.2 g/gDCW for the simulations, and higher for the actual data (includes all additional proteins quantified in experiment, but not carrying flux in model simulations).

df_prot_comp %>% group_by(condition, replicate) %>% 
  filter(predicted_mass_g_gDCW != 0) %>%
  summarize(sum(mass_g_gDCW), sum(predicted_mass_g_gDCW))

`summarise()` has grouped output by 'condition'. You can override using the `.groups` argument.

---

The total predicted protein mass is lower than the measured protein mass. Therefore the following section quantifies discrepancies between model predicted and actually measured abundances. First we can inspect the top N reactions with highest average predicted protein abundance. The ratio of predicted divided by measured mass indicates that a handful of proteins are predicted to be more than 10 fold abundant compared to the measured abundance. This points towards fluxes being erroneously predicted too high for particular reactions, or k_app values being estimated too low for the estimated flux.

However, the largest discrepancies arise from under-estimation of proteins, the main cause being that the model predicts the optimal abundance for each enzyme to carry a certain flux. If fluxes are drastically reduced due to strong substrate limitation, the minimal required protein abundance to optimize growth will be much lower than the measured abundance. A bacterial cell on the other hand can not fully reduce its proteome but instead ‘suspends’ inactive enzymes.

6.2 Change of machinery proteins with growth rate

The simulation and measurement data was prepared and merged by condition in the previous sections. Now it can be plotted to e.g. compare protein allocation over growth rate. Interestingly, we see that model predictions are quite accurately reflecting the range of protein allocation for the four different machineries, see following paragraphs. This is a good confirmation of the model’s predictive power, given that the rates of these machineries were not fitted from data but taken purely from literature. There are however some deviations from the predicted ‘optimal’ proteome:

• the most important machine is the ribosome. Prediction and experiment show a very similar increase of ribosome abundance with growth rate (slope of linear model), but the intersection (amount of unused ribosomes at zero growth) is much higher in experiment
• chaperones show an inverse proportional relationship with growth rate contrary to model prediction. Do (some of?) these proteins have another role than just folding, like stress response?
• transcription sector is quite stable, however, predicted enzyme mass is much lower than measured (underestimated by 1 order of magnitude , 5x10^-3 vs 5x10^-4 g/gDCW)
• replication sector is heavily underestimated by 3 orders of magnitude (10^-3 vs 10^-6 g/gDCW)

The protein allocation for the two low-abundant machines, replication and transcription is very hard to see. We can update the Y-axis limits for this plot and re-plot it with 10-fold ‘zoom’. Abundance of replication machinery increases with growth rate from 0.0010 to a maximum of 0.0015 g/gDCW, an increase by 20-50%. For transcription, abundance increases with growth rate from 0.006 to 0.008 g/gDCW, an increase of 33%.

df_machinery_util <- df_machinery %>%
  
  # remove 1 non-quantified replicate
  filter(!(substrate == "fructose" & growth_rate == 0.1 & replicate == "R2")) %>%
  
  # summing up protein mass over all conditions
  group_by(machine, substrate, growth_rate, replicate) %>%
  summarize(
    `mass [g/gDCW]` = sum(mass_g_gDCW, na.rm = TRUE),
    `predicted mass [g/gDCW]` = sum(predicted_mass_g_gDCW, na.rm = TRUE),
    utilization = `predicted mass [g/gDCW]`/`mass [g/gDCW]`
  ) %>% ungroup %>%
  
  # replace 0 with NA and reorder factors
  mutate(across(matches("mass"), ~ na_if(.x, 0))) %>%
  mutate(machine = machine %>% factor(., unique(.)[c(2,4,3,1)]))

`summarise()` has grouped output by 'machine', 'substrate', 'growth_rate'. You can override using the `.groups` argument.

# plot mass and utilization
plot_machinery <-
  
  xyplot(`mass [g/gDCW]` + `predicted mass [g/gDCW]` ~ factor(growth_rate) | 
    machine * substrate, df_machinery_util,
    par.settings = custom.colorblind(), auto.key = list(columns = 2, cex = 0.7),
    xlab = expression("µ [h"^-1*"]"),
    ylab = expression("m"[protein]*" [g gDCW"^-1*"]"),
    scales = list(alternating = FALSE, x = list(at = c(2,4))),
    as.table = TRUE, ylim = c(-0.01, 0.12),
    between = list(x = 0.5, y = 0.5), pch = 19, lwd = 2,
    panel = function(x, y, ...) {
      panel.grid(h = -1, v = -1, col = grey(0.9))
      panel.errbars(x, y, ewidth = 0, ...)
      panel.superpose(x, y, ...)
    }, panel.groups = function(x, y, ...) {
      panel.lmline(x, y, ...)
    }
  )

6.3 Under-utilization of machinery proteins

Here, we determine the underutilized machinery fraction by taking the ratio of simulated optimal enzyme abundance and experimentally measured abundance, over all four central dogma machines.

plot_machinery_util <- df_machinery %>%
  
  # remove 1 non-quantified replicate
  filter(!(substrate == "fructose" & growth_rate == 0.1 & replicate == "R2")) %>%
  
  # summarizing protein utilization for all machines
  mutate(machine = "machines") %>%
  group_by(machine, substrate, growth_rate, replicate) %>%
  summarize(
    mass = sum(mass_g_gDCW, na.rm = TRUE),
    predicted_mass = sum(predicted_mass_g_gDCW, na.rm = TRUE),
    utilization = predicted_mass/(mass)
  ) %>%
  
  # plot. use alternating = 2 to switch axis to right side
  xyplot(utilization*100 ~ factor(growth_rate) | machine * substrate, .,
    par.settings = custom.colorblind(), 
    scales = list(x = list(at = c(2,4)), y = list(alternating = FALSE)),
    xlab = expression("µ [h"^-1*"]"), ylab = "",
    key = simpleKey("% utilization", cex = 0.7),
    as.table = TRUE, between = list(x = 0.5, y = 0.5),
    pch = 19, lwd = 2, layout = c(1, 4), ylim = c(-10, 110),
    panel = function(x, y, ...) {
      panel.grid(h = -1, v = -1, col = grey(0.9))
      x_mean = unique(x); y_mean = tapply(y, x, mean)
      panel.xyarea(c(0, x_mean, 6), c(0, y_mean, tail(y_mean, 1)), 
        lty = 0, col = grey(0.6, alpha = 0.5), ...)
      panel.errbars(x, y, ewidth = 0, col = grey(0.5), ...)
    }
  ) %>% useOuterStrips

`summarise()` has grouped output by 'machine', 'substrate', 'growth_rate'. You can override using the `.groups` argument.

print(useOuterStrips(plot_machinery), position = c(0, 0, 0.77, 1.027), more = TRUE)
print(plot_machinery_util, position = c(0.71, 0, 1.02, 1.027))

grid::grid.text(c("B", "C"), x = c(0.02, 0.75), y = c(0.97,0.97))

7 Under-utilization of enzymes

7.1 Global trends in utilization

We can subdivide enzymes by utilization and plot average mass, average variability, and average essentiality. Average essentiality here means that every reaction that is associated with at least one probably essential genes is probably essential, and every reactions with at least one essential gene is counted as essential (taking precedence over ‘probably essential’).

df_util <- df_prot_comp %>%
  
  # filter reactions without protein measurement or substrate
  filter(!(is.na(mass_g_gDCW) | is.na(substrate)), growth_rate == 0.25) %>%
  
  # first calculate utilization per condition
  group_by(model_group_basic, reaction_id, reaction_name, substrate) %>%
  summarize(.groups = "drop_last",
    mass_g_gDCW = mean(mass_g_gDCW, na.rm = TRUE),
    predicted_mass_g_gDCW = mean(predicted_mass_g_gDCW, na.rm = TRUE),
    utilization = predicted_mass_g_gDCW/mass_g_gDCW
  ) %>%
  
  # determine average mass, sd and utilization per reaction
  summarize(.groups = "drop_last",
    mean_mass = mean(mass_g_gDCW, na.rm = TRUE),
    sd_mass = sd(mass_g_gDCW, na.rm = TRUE),
    cv_mass = sd_mass/mean_mass,
    utilization = mean(utilization, na.rm = TRUE)
  ) %>%
  
  # manual 'clustering' based on protein vs growth rate slope
  mutate(utilization_group = case_when(
    utilization <= 0.33 ~ "low",
    between(utilization, 0.33, 0.66) ~ "moderate",
    utilization > 0.66 ~ "high",
  ) %>% factor(., c("low", "moderate", "high")))


# import essentiality from TnSeq/BarSeq data
df_util <- read_csv("../data/output/essentiality_escher.csv", col_types = cols()) %>%
  mutate(essentiality = essentiality %>% recode(
    `0` = "not essential", `1` = "probably essential", `2` = "essential")) %>%
  mutate(essentiality = case_when(
    locus_tag %in% read_csv("../data/output/essentiality_fructose.csv", col_types = cols())$locus_tag ~ "essential",
    TRUE ~ essentiality
  )) %>%
  # merge with utilization
  inner_join(df_model_reactions %>% rename(locus_tag = genes)) %>%
    group_by(reaction_id) %>% summarize(essentiality = max(essentiality, na.rm = TRUE) %>%
    recode(`0` = "not essential", `1` = "probably essential", `2` = "essential")) %>%
  right_join(df_util)

Joining, by = "locus_tag"
Joining, by = "reaction_id"

Test if mean enzyme mass and enzyme variability (CV) are normal-distributed. This is important for choosing the statistical test for significance. We use the Shapiro-Wilks test for normality.

df_util %>% group_by(utilization_group) %>%
  summarize(
    Shapiro_Wilk_pvalue_mass = shapiro.test(log10(mean_mass))$p.value,
    Shapiro_Wilk_pvalue_CV = shapiro.test(log10(cv_mass))$p.value
  )

Now plot average mass in g/gDCW, CV of mass in g/gDCW and essentiality (score) from TnSeq experiments. As some of the data is not normally distributed, we apply a non-parametric significance test instead - here, for comparison of two samples, the best alternative to Student’s t-test is the Mann-Whitney U-test. For one sample comparisons, the non-parametric alternative to t-test would be the Wilcoxon Rank Sum test.

boxplot_colorblind <- custom.colorblind()
boxplot_colorblind$box.rectangle$lwd = 1.5
boxplot_colorblind$box.umbrella$lwd = 1.5


plot_underutil_mass <- xyplot(log10(mean_mass) ~ utilization_group,
    df_util,
    par.settings = boxplot_colorblind,
    horizontal = FALSE, do.out = FALSE, ylim = c(-8, -1),
    xlab = "utilization", ylab = expression("m"[protein]*" [log"[10]*" g gDCW"^-1*"]"),
    scales = list(x = list(cex = 0.7)),
    panel = function(x, y, ...) {
      panel.grid(h = -1, v = -1, col = grey(0.9))
      panel.violin(x, y, lwd = 0, col = grey(0.6, 0.5), box.width = 0.7, ...)
      panel.bwplot(x, y, pch = "|", box.width = 0.3, ...)
      panel.pvalue(x, y, cex = 0.7, fixed_pos = -0.9, verbose = TRUE, method = "wilcox.test", ...)
      panel.text(unique(x), -7.5, labels = paste0("n=", table(x)), cex = 0.6)
    }
  )

plot_underutil_cv <- xyplot(log10(cv_mass) ~ utilization_group,
    df_util,
    par.settings = boxplot_colorblind,
    horizontal = FALSE, do.out = FALSE, ylim = c(-1.6,0.6),
    xlab = "utilization", ylab = expression("CV m"[protein]*" [log"[10]*"]"),
    scales = list(x = list(cex = 0.7)),
    panel = function(x, y, ...) {
      panel.grid(h = -1, v = -1, col = grey(0.9))
      panel.violin(x, y, lwd = 0, col = grey(0.6, 0.5), box.width = 0.7, ...)
      panel.bwplot(x, y, pch = "|", box.width = 0.3, ...)
      panel.pvalue(x, y, cex = 0.7, fixed_pos = 0.65, verbose = TRUE, method = "wilcox.test", ...)
      panel.text(unique(x), -1.45, labels = paste0("n=", table(x)), cex = 0.6)
    }
  )

plot_underutil_es <- xyplot(n_reactions ~ utilization_group,
    df_util %>% group_by(utilization_group, essentiality) %>% summarize(n_reactions = n()),
    par.settings = boxplot_colorblind,
    col = custom.colorblind()$superpose.polygon$col[c(5,4,2)],
    groups = factor(essentiality, c("not essential", "probably essential", "essential")),
    xlab = "utilization", ylab = "number of reactions", ylim = c(-20, 620),
    ewidth = 0.1, lwd = 2, scales = list(x = list(cex = 0.7)),
    panel = function(x, y, ...) {
      panel.grid(h = -1, v = -1, col = grey(0.9))
      panel.barplot(x, y, beside = TRUE, ...)
      panel.key(..., points = FALSE, corner = c(0.95, 0.95), cex = 0.7)
    }
  )

`summarise()` has grouped output by 'utilization_group'. You can override using the `.groups` argument.

print(plot_underutil_mass, position = c(0,0,0.35,1.05), more = TRUE)

print(plot_underutil_cv, position = c(0.32,0,0.68,1.05), more = TRUE)

print(plot_underutil_es, position = c(0.64,0,1,1.05))
grid::grid.text(c("A", "B", "C"), x = c(0.02,0.35,0.65), y = c(0.94,0.94,0.94))

We can try to combine several types of information for each reaction in concise mini-figures:

• the essentiality (simply by color around reaction name)
• the average protein mass associated with the reaction
• the change over conditions (CV)
• the utilization/saturation

# enzymes grouped by pathway
list_ccm_enzymes <- list(
  cbb = c("PGK", "TPI", "GAPD", "FBA", "FBP", "TKT1", "TKT2", "TAh", "RPE", "RPI", "PRUK", "RBPC", "FDH"),
  ed = c("PGI", "G6PDH2r", "PGL", "EDD", "EDA"),
  pyr = c("PGM", "ENO", "PYK", "PDH1", "PDH2", "PDH3", "PC", "PPC", "ME1", "ME2", "CS"),
  tca = c("ACONT1", "ACONT2", "ICDHx", "ICDHyr", "AKGDH", "SUCOAS", "SUCDi", "SUCD1", "FUM", "MDH", "MALS")
)

plot_minifigs <- df_prot_comp %>%
  filter(reaction_id %in% unlist(list_ccm_enzymes),
    growth_rate == 0.25) %>%
  group_by(reaction_id, substrate) %>%
  # calculate utilization
  mutate(utilization = predicted_mass_g_gDCW/mass_g_gDCW) %>%
  summarize(utilization = mean(utilization), mean_mass = mean(mass_g_gDCW)) %>%
  # rescale all to new range between 0 and 1
  ungroup %>%
  mutate(across(matches("mass|utilization$"), scales::rescale)) %>%
  pivot_longer(all_of(c("mean_mass", "utilization"))) %>%
  mutate(name = factor(name, c("mean_mass", "utilization"))) %>%
  
  xyplot(value ~ name | reaction_id, .,
    groups = substrate, par.settings = custom.colorblind(),
    xlab = "", ylab = "", fill_alpha = 1,
    ewidth = 0.09, lty = 1.5, as.table = TRUE,
    scales = list(draw = FALSE), ylim = c(-0.1, 1.0),
    between = list(x = 0.5, y = 0.5),
    panel = function(x, y, ...) {
      panel.rect(0.5,-0.1,2.5,1.1, col = "white", border = "transparent")
      panel.rect(0.5,-0.1,1.5,1.0, col = grey(0.9), border = "transparent")
      panel.barplot(x, y, beside = TRUE, ...)
    }
  )

`summarise()` has grouped output by 'reaction_id'. You can override using the `.groups` argument.

print(plot_minifigs)

7.2 Central carbon metabolism, detailed utilization per enzyme

Analogously to the above analysis, we can plot protein abundance and utilization for all enzymes of central carbon metabolism (from Figure 3 D).

plot_enz_ccm <- lapply(list_ccm_enzymes, function(lst) {
  df_prot_comp %>%
  filter(reaction_id %in% lst, !is.na(condition)) %>%
  mutate(reaction_id = factor(reaction_id, lst)) %>%
  group_by(reaction_id) %>% mutate(mass_g_gDCW = scales::rescale(mass_g_gDCW)) %>%
  
  xyplot(mass_g_gDCW ~ factor(growth_rate) | reaction_id, .,
    par.settings = custom.colorblind(),
    groups = substrate, xlab = "", #expression("µ [h"^-1*"]"),
    #ylab = expression("m"[protein]*" [g gDCW"^-1*"]"),
    ylab = "relative abundance",
    layout = c(8, ifelse(length(lst) > 8, 2, 1)),
    ewidth = 0, lwd = 2, as.table = TRUE,
    between = list(x = 0.5, y = 0.5),
    scales = list(alternating = FALSE, x = list(at = c(2, 4))),
    panel = function(x, y, ...) {
      panel.grid(h = -1, v = -1, col = grey(0.9))
      panel.errbars(x, y, ...)
    }
  )
})

print(plot_enz_ccm[[1]], position = c(0.0,0.67,1,1), more = TRUE)
print(plot_enz_ccm[[2]], position = c(0,0.51,1,0.75), more = TRUE)

print(plot_enz_ccm[[3]], position = c(0,0.26,1,0.59), more = TRUE)
print(plot_enz_ccm[[4]], position = c(0,0,1,0.34))

grid::grid.text(c("CBB", "ED", "PYR", "TCA"), x = c(0.03, 0.03, 0.03, 0.03), y = c(0.97, 0.71, 0.53, 0.3))

And we can add utilization for the different enzymes and conditions too.

plot_util_ccm <- lapply(list_ccm_enzymes, function(lst) {
  df_prot_comp %>%
  filter(reaction_id %in% lst, !is.na(condition)) %>%
  mutate(reaction_id = factor(reaction_id, lst)) %>%
  # calculate utilization
  mutate(percent_utilization = predicted_mass_g_gDCW/mass_g_gDCW*100) %>%
  
  xyplot(percent_utilization ~ factor(growth_rate) | reaction_id, .,
    par.settings = custom.colorblind(),
    groups = substrate,
    xlab = "", ylab = "% utilization", ylim = c(-20, 200),
    layout = c(8, ifelse(length(lst) > 8, 2, 1)),
    ewidth = 0, lwd = 2, as.table = TRUE,
    between = list(x = 0.5, y = 0.5),
    scales = list(alternating = FALSE, x = list(at = c(2, 4))),
    panel = function(x, y, ...) {
      panel.grid(h = -1, v = -1, col = grey(0.9))
      panel.errbars(x, y, ...)
    }
  )
})

print(plot_util_ccm[[1]], position = c(0.0,0.67,1,1), more = TRUE)
print(plot_util_ccm[[2]], position = c(0,0.51,1,0.75), more = TRUE)

print(plot_util_ccm[[3]], position = c(0,0.26,1,0.59), more = TRUE)
print(plot_util_ccm[[4]], position = c(0,0,1,0.34))

grid::grid.text(c("CBB", "ED", "PYR", "TCA"), x = c(0.03, 0.03, 0.03, 0.03), y = c(0.97, 0.71, 0.53, 0.3))

Finally we can take a look on trends with growth rate, by fitting a linear model for each enzyme and condition over growth rate, and plotting the R or R squared per condition. The R (correlation coefficient) should indicate the trend, i.e. increasing or decreasing with growth rate for a certain substrate.

plot_lm_ccm <- lapply(list_ccm_enzymes, function(lst) {
  
  df_prot_comp %>%
  filter(reaction_id %in% lst, !is.na(condition)) %>%
  mutate(reaction_id = factor(reaction_id, lst)) %>%
  group_by(reaction_id) %>% mutate(mass_g_gDCW = scales::rescale(mass_g_gDCW)) %>%
  
  # fit linear model
  group_by(reaction_id, substrate) %>%
  summarize(
    lin_reg_slope = summary(lm(mass_g_gDCW ~ growth_rate))$coef[2],
    lin_reg_pvalue = summary(lm(mass_g_gDCW ~ growth_rate))$coef[8]
  ) %>%
  mutate(lin_reg_slope = lin_reg_slope %>% 
    if_else(. > 2.5, 2.5, .) %>% if_else(. < -2.5, -2.5, .)) %>%
  
  levelplot(lin_reg_slope ~ reaction_id * factor(substrate, unique(substrate)[c(3,2,1,4)]), .,
    par.settings = custom.colorblind(), at = seq(-2.5, 2.5, by = 0.25),
    col.regions = colorspace::diverge_hcl(20), ylab = "",
    as.table = TRUE, xlab = "",
    scales = list(x = list(cex = 0.65))
  )
})

`summarise()` has grouped output by 'reaction_id'. You can override using the `.groups` argument.
`summarise()` has grouped output by 'reaction_id'. You can override using the `.groups` argument.
`summarise()` has grouped output by 'reaction_id'. You can override using the `.groups` argument.
`summarise()` has grouped output by 'reaction_id'. You can override using the `.groups` argument.

print(plot_lm_ccm[[1]], position = c(0.0,0.68,1,1), more = TRUE)
print(plot_lm_ccm[[2]], position = c(0,0.46,1,0.77), more = TRUE)

print(plot_lm_ccm[[3]], position = c(0,0.24,1,0.55), more = TRUE)
print(plot_lm_ccm[[4]], position = c(0,0,1,0.32))

grid::grid.text(c("CBB", "ED", "PYR", "TCA"), x = c(0.03, 0.03, 0.03, 0.03), y = c(0.97, 0.73, 0.5, 0.28))

7.3 Example of underutilization of enzymes from Cbb operon

Similar to machinery proteins we can also follow the simulation and actual protein abundance of all (detected) enzymes. We can see that some protein abundances nicely follow a growth-rate dependent manner, in line with predictions of higher flux. One such example are Calvin cycle enzymes for formatotrophic growth. However, for the same enzymes we see that no abundance is predicted under heterotrophic conditions because of missing flux through the pathway, but in fact the proteins are expressed in high abundance, often in a growth-rate dependent manner.

cbb_selected <- c("PGK", "GAPD", "FBA", "FBP", "TKT1", "PRUK", "RBPC")
# "TKT2", "TAh", "RPE", "RPI"

# summary table of selected Calvin cycle genes
df_calvin <- df_prot_comp %>%
  
  # remove 1 non-quantified replicate
  filter(!(substrate == "fructose" & growth_rate == 0.1 & replicate == "R2")) %>%
  
  # select only subset of reactions
  filter(reaction_id %in% cbb_selected) %>%
  mutate(reaction_id = factor(reaction_id, cbb_selected)) %>%
  rename(`mass [g/gDCW]` = mass_g_gDCW, `predicted mass [g/gDCW]` = predicted_mass_g_gDCW) %>%
  mutate(utilization = `predicted mass [g/gDCW]`/`mass [g/gDCW]`)


# plot mass and utilization
plot_calvin <- #doubleYScale(use.style = FALSE, under = TRUE,
  
  xyplot(`mass [g/gDCW]` + `predicted mass [g/gDCW]` ~ factor(growth_rate) | 
    reaction_id * substrate, df_calvin,
    par.settings = custom.colorblind(), auto.key = list(columns = 2, cex = 0.7),
    xlab = expression("µ [h"^-1*"]"), ylim = c(-0.003, 0.033), #c(-0.002, 0.032),
    ylab = expression("m"[protein]*" [g gDCW"^-1*"]"),
    scales = list(alternating = FALSE, x = list(at = c(2,4)),
      y = list(at = c(0, 0.01, 0.02, 0.03))), as.table = TRUE,
    between = list(x = 0.5, y = 0.5), lwd = 2, 
    panel = function(x, y, ...) {
      panel.grid(h = -1, v = -1, col = grey(0.9))
      panel.errbars(x, y, ewidth = 0, ...)
      panel.superpose(x, y, ...)
    }, panel.groups = function(x, y, ...) {
      panel.lmline(x, y, ...)
    }
  )

print(useOuterStrips(plot_calvin))

It is clear from the previous analysis that the cells maintain proteins even if they are not or only marginally utilized. The actually utilized proteome is minimal under strong substrate limitation. The next section will try to quantify the under-utilization of enzymes by comparing the minimal protein requirement (simulation) versus experimentally determined protein abundance. It is sufficient to determine the underutilized protein fraction (of all utilized proteins).

plot_calvin_util <- df_prot_comp %>%
  
  # filter set of enzymes
  filter(reaction_id %in% c("PGK", "GAPD", "FBA", "FBP", "TKT1", "PRUK", "RBPC")) %>%
  
  # remove 1 non-quantified replicate
  filter(!(substrate == "fructose" & growth_rate == 0.1 & replicate == "R2")) %>%
  
  # calculate protein utilization
  group_by(substrate, growth_rate, replicate) %>%
  summarize(
    mass = sum(mass_g_gDCW, na.rm = TRUE),
    predicted_mass = sum(predicted_mass_g_gDCW, na.rm = TRUE),
    percent_utilization = predicted_mass/(mass)
  ) %>%
  
  # add pseudo enzyme group
  mutate(enzyme = "CBB cycle") %>%
  
  # optionally show some summary statistics
  #summarize(across(matches("mass|util"), mean)) %>%
  #mutate(unutil_mass = mass-predicted_mass) %>%
  #summarize(max(unutil_mass)/0.68*100) #for percent: /0.68*100
  
  # plot
  xyplot(percent_utilization*100 ~ factor(growth_rate) | enzyme * substrate, .,
    par.settings = custom.colorblind(), ylim = c(-10, 110),
    scales = list(alternating = FALSE, x = list(at = c(2,4))),
    xlab = expression("µ [h"^-1*"]"), ylab = "",
    as.table = TRUE, between = list(x = 0.5, y = 0.5),
    lwd = 2, key = simpleKey("% utilization", cex = 0.7),
    panel = function(x, y, ...) {
      panel.grid(h = -1, v = -1, col = grey(0.9))
      x_mean = unique(x); y_mean = tapply(y, x, mean)
      panel.xyarea(c(0, x_mean, 6), c(0, y_mean, tail(y_mean, 1)), 
        lty = 0, col = grey(0.6, alpha = 0.5), ...)
      panel.errbars(x, y, ewidth = 0, col = grey(0.5), ...)
    }
  ) %>% useOuterStrips

`summarise()` has grouped output by 'substrate', 'growth_rate'. You can override using the `.groups` argument.

print(plot_calvin_util)

---

A final control is to check which copies of the same Cbb cycle enzymes are more abundant, or undergo stronger changes. For example, PRUK and Rubisco (RBPC) are only present on the two copies of the Calvin cycle operon, on pHG1 and Chromosome 2. But other glycolysis related enzymes have one canonical copy on chromosome 1, which seems to encode mostly “housekeeping” functions. We can look at protein abundance of the three single gene loci for each enzyme (chromosome 1, 2, and pHG1 megaplasmid). It’s important to keep in mind that most peptides for cbb genes can not be distinguished between chromosome 2 and pHG1, therefore the similar expression pattern.

plot_calvin_loci <- df_model_reactions %>% filter(reaction_id %in% c("PGK", "GAPD", "FBA", "FBP", "TKT1")) %>%
  select(reaction_id, genes) %>% rename(locus_tag = genes) %>%
  filter(locus_tag != "H16_B0278") %>%
  inner_join(Ralstonia_eutropha) %>%
  mutate(chromosome = case_when(
    grepl("H16_A", locus_tag) ~ "chromosome 1",
    grepl("H16_B", locus_tag) ~ "chromosome 2",
    grepl("PHG", locus_tag) ~ "pHG1",
  )) %>%
  
  xyplot(mass_g_gDCW ~ factor(growth_rate) | reaction_id * chromosome, .,
    groups = substrate, layout = c(5, 3),
    par.settings = custom.colorblind(), lwd = 2,
    scales = list(alternating = FALSE),
    xlab = expression("µ [h"^-1*"]"),
    ylab = expression("m"[protein]*" [g gDCW"^-1*"]"),
    as.table = TRUE, between = list(x = 0.5, y = 0.5),
    panel = function(x, y, ...) {
      panel.grid(h = -1, v = -1, col = grey(0.9))
      panel.errbars(x, y, ewidth = 0, type = c("p", "l"), ...)
      panel.key(..., cex = 0.7)
    }
  )

Joining, by = "locus_tag"

print(plot_calvin_loci)

8 Expression and utilization of PHB synthesis enzymes

The pathway for PHB synthesis consists mainly of three enzymes, each one with multiple genes annotated. The first is Acetyl-CoA C-acetyltransferase (phaA), making acetoacetyl-CoA, the second one is Acetoacetyl-CoA reductase (phaB), making 3-hydroxy-butyryl-CoA, and the third is the PHB-synthase condensing 3HB monomers to PHB polymer (phaC). The mass for each enzyme is the sum of the individual protein masses allocated to this reaction.

plot_phb_mass <- df_prot_comp %>%
  filter(reaction_id %in% c("PHAS", "ACACT1r", "AACOAR")) %>%
  mutate(reaction_id = recode(reaction_id,
    PHAS = "PHAS (phaC)", ACACT1r = "ACACT1r (phaA)", AACOAR = "AACOAR (phaB)")) %>%
  
  xyplot(mass_g_gDCW ~ factor(growth_rate) | 
      factor(reaction_id, unique(reaction_id)[c(2,1,3)]), .,
    par.settings = custom.colorblind(),
    groups = substrate, xlab = "",
    ylab = expression("m"[protein]*" [g gDCW"^-1*"]"),
    layout = c(3, 1), ewidth = 0, lwd = 2, as.table = TRUE,
    between = list(x = 0.5, y = 0.5), auto.key = list(columns = 2),
    scales = list(alternating = FALSE, x = list(at = c(2, 4))),
    panel = function(x, y, ...) {
      panel.grid(h = -1, v = -1, col = grey(0.9))
      panel.errbars(x, y, ...)
    }
  )

plot_phb_util <- df_prot_comp %>%
  filter(reaction_id %in% c("PHAS", "ACACT1r", "AACOAR")) %>%
  mutate(percent_utilization = predicted_mass_g_gDCW/mass_g_gDCW*100) %>%
  mutate(reaction_id = recode(reaction_id, 
    PHAS = "PHAS (phaC)", ACACT1r = "ACACT1r (phaA)", AACOAR = "AACOAR (phaB)")) %>%
  
  xyplot(percent_utilization ~ factor(growth_rate) | 
      factor(reaction_id, unique(reaction_id)[c(2,1,3)]), .,
    par.settings = custom.colorblind(),
    groups = substrate,
    xlab = expression("µ [h"^-1*"]"),
    ylab = "% utilization", ylim = c(-20, 320),
    layout = c(3, 1),
    ewidth = 0, lwd = 2, as.table = TRUE,
    between = list(x = 0.5, y = 0.5),
    scales = list(alternating = FALSE, x = list(at = c(2, 4))),
    panel = function(x, y, ...) {
      panel.grid(h = -1, v = -1, col = grey(0.9))
      panel.errbars(x, y, ...)
    }
  )

print(plot_phb_mass, position = c(0,0.42,1,1), more = TRUE)
print(plot_phb_util, position = c(0.04,0,1,0.55))

9 Yield-growth rate tradeoff when Cupriavidus re-assimilates CO2?

It was experimentally observed before that R. eutropha expresses Rubisco and other cbb-operon located genes even on growth on fructose or other heterotrophic substrates. Bowien et al., 1990, show that Rubisco activity can be found on growth on fructose, but not on pyruvate. Dangel & Tabita, 2015, review the regulation by CbbR type regulators among others also in R. eutropha and mention that citrate also leads to activation of cbb expression. They hypothesize that ribulose bisphosphate (RuBP) is an activating effector while other central metabolism intermediates such as phosphoenolpyruvate (PEP) are repressing effector molecules. It was speculated in Bowien et al. that regulation by cbbR might actually be a repression-derepression rather than activation, which means that the default state would be a bound cbbR repressing the cbb operon. However it became clear that this is not the case. Shimizu et al., 2015, knocked the cbbR gene out and the result was reduced expression of cbb genes by 100 fold. This proves that cbbR is a required activator and not a repressor of cbb, otherwise cbbR deletion would lead to constitutive activation of cbb expression.

The same group speculated in their study that the additional Rubisco expression could have a benefit for carbon yield, specifically for product yield of PHB. They show that PHB from the WT contains slightly more 13C labeled mass (and total mass) than the cbbR and cbbS/L knockouts. This means that the cell would have a (product or biomass) yield advantage by Rubisco expression on fructose.

The following section tests the hypothesis of a yield-growth rate tradeoff with the resource allocation model. First simulation data from the RBA model is imported

#adjust path
mixotroph_dir <- gsub("substrate_limitation", "mixotrophy", simulation_dir)

# read simulation results
df_mixflux <- read_rba_result(list.files(mixotroph_dir, pattern = "fluxes_.*.tsv$", full.names = TRUE))
df_mixmacr <- read_rba_result(list.files(mixotroph_dir, pattern = "macroprocesses_.*.tsv", full.names = TRUE))

# rename column
df_mixflux <- df_mixflux %>% rename(CO2_refixation = sim_run) %>% filter(CO2_refixation <= 5)
df_mixmacr <- df_mixmacr %>% rename(CO2_refixation = sim_run) %>% filter(CO2_refixation <= 5)

The next task is to compare yield, growth rate, CO₂ emission and other fluxes for a range of simulations where Rubisco was forced to re-fix emitted CO₂ from growth on fructose. Simulations were performed for increasing flux through Rubisco from 0 to 5 mmol gDCW^-1 h^-1.

plot_mixo_mu <- lapply(c("mu", "yield"), function(keys) {
  
  if (keys == "mu") {
    ylabel <- expression("µ [h"^-1*"]")
    ylim <- c(0.001, 0.3)
  } else {
    ylabel <- expression("Y [gDCW gS"^-1*"]")
    ylim <- c(0.001, 0.399)
  }
  
  
  xyplot(value ~ factor(CO2_refixation),
    filter(df_mixmacr, key == keys, carbon_conc == 1.25),
    type = "b", lwd = 2, par.settings = custom.colorblind(), 
    xlab = expression("CO"[2]*" uptake [mmol h"^-1*" gDCW"^-1*"]"),
    ylab = ylabel, ylim = ylim,
    panel = function(x, y, z, ...) {
      panel.grid(h = -1, v = -1, col = grey(0.9))
      panel.xyplot(x, y, ...)
    }
  )
})

print(plot_mixo_mu[[1]], position = c(0,0,0.53,1), more = TRUE)
print(plot_mixo_mu[[2]], position = c(0.47,0,1,1))

There is no increase in growth rate or yield with additional CO₂ fixation according to the RBA model. The yield is in fact calculated based on the growth rate µ and the substrate uptake rate qS. Yield and growth rate depend on each other in the relation Y [gDCW/gS] = µ [h-1] / qS [gS gDCW-1 h-1].

Since only the initial substrate concentration is constrained, the model could predict a lower substrate uptake rate or lower CO₂ emission per biomass in order to reach a yield increase. However this was not the case under any simulation, see below for details of specific metabolites.

The model seems to predict a high-yield phenotype already. If CO₂ re-fixation would have had an advantage for growth, it might have been detected earlier depending on the protein costs for the respective pathway. Before a final verdict, we can follow the fate of the fixed CO₂ through the metabolism.

plot_mixo_enz <- xyplot(abs(value) ~ as.factor(CO2_refixation) | 
    factor(key, unique(key)), layout = c(2, 2),
  filter(df_mixflux, key %in% c("CO2t", "EDD", "CS", "O2t"), carbon_conc == 1.25),
  par.settings = custom.colorblind(),
  type = "b", lwd = 2, between = list(x = 0.5, y = 0.5), as.table = TRUE,
  xlab = expression("CO"[2]*" uptake [mmol h"^-1*" gDCW"^-1*"]"),
  ylab = expression("flux [mmol h"^-1*"gDCW"^-1*"]"),
  scales = list(alternating = FALSE),
  panel = function(x, y, z, ...) {
    panel.grid(h = -1, v = -1, col = grey(0.9))
    panel.xyplot(x, y, ...)
  }
)

print(plot_mixo_enz)

9.1 Conclusion

It becomes clear that yield and growth rate decrease, and not increase, with additional CO2 fixation, because:

• energy requirement for CO2 fixation leads to lower flux through ED pathway, but higher flux through TCA
• this is in order to generate the required NADH/ATP for CO2 fixation
• respiration and O2 consumption also increases with forced mixotrophic growth
• there is no net reduction of CO2 emission. In fact cells emit more CO2 even when they fix some of it due to increased energy requirement
• storing some of the fixed CO2 as PHB would not increase biomass yield as additional energy requirement still persists
• cells should not be able to make more PHB using this strategy, regardless of biomass yield. If acetyl-CoA is drained for PHB, even less energy is made available through TCA and OXPHOS.

10 Miniature RBA flux maps

Install my small R package fluctuator to overlay fluxes on SVG metabolic maps. The template map is a simplified metabolic map of C. necator’s central metabolism.

First we plot a map of the mixotrophic simulation. Then we use the same work-flow and template to plot metabolic flux maps for growth on the three main substrates, fructose, succinate and formate (maximum tested growth rate, 0.25 h^-1).

library(fluctuator)

# import map
SVG_map <- read_svg("../data/simulation/central_metabolism.svg")
#get_attributes(SVG_map, node = "key_01") %>% pull(style)

# prepare data: mixotrophic growth
df_mixflux <- df_mixflux %>% filter(CO2_refixation == 3) %>%
  # add fluxes for legend entries
  add_row(key = c("key_01", "key_12", "key_24", "key_48"), value = c(1,2,4,8)) %>%
  # add stroke width to simulation
  mutate(stroke_width = 0.2 + 0.2*sqrt(abs(value))) %>%
  inner_join(select(SVG_map@summary, label), by = c("key" = "label")) %>%
  distinct

# prepare data: carbon-limited growth
df_carbflux <- df_flux %>%
  filter(simulation %in% c(
    "for_0.487_nh4_18.7_10", "fru_1.223_nh4_18.7_0", "succ_0.662_nh4_18.7_5")) %>%
  # add fluxes for legend entries
  add_row(key = rep(c("key_01", "key_12", "key_24", "key_48"), each = 3),
    value = rep(c(1,2,4,8), each = 3), carbon_source = rep(c("for", "fru", "succ"), 4)) %>%
  # add stroke width to simulation
  mutate(stroke_width = 0.2 + 0.2*sqrt(abs(value))) %>%
  inner_join(select(SVG_map@summary, label), by = c("key" = "label")) %>%
  distinct

Plot all four conditions in a loop.

lapply(
  list(
    mutate(df_mixflux, carbon_source = "mix"),
    filter(df_carbflux, carbon_source == "fru"),
    filter(df_carbflux, carbon_source == "for"),
    filter(df_carbflux, carbon_source == "succ")
  ), FUN = function(df) {
    
    # reload base map for every instance, as the SVG is altered
    # even without assignment (python style, untypical in R)
    SVG_map <- read_svg("../data/simulation/central_metabolism.svg")
    
    # modify map
    SVG2 <- set_attributes(SVG_map,
      node = df$key, attr = "style",
      pattern = "stroke-width:[0-9]+\\.[0-9]+",
      replacement = paste0("stroke-width:", df$stroke_width))
    
    # turn non-zero fluxes into darker color
    SVG2 <- set_attributes(SVG2,
      node = filter(df, value != 0) %>% pull(key), attr = "style",
      pattern = "stroke:#b3b3b3",
      replacement = "stroke:#808080")
    
    # rescale arrow heads
    SVG2 <- set_attributes(SVG2, 
      node = grep("marker", SVG2@summary$id, value = TRUE),
      node_attr = "id",
      attr = "transform",
      pattern = "scale\\(0.2\\)",
      replacement = "scale(0.15)")
    
    SVG2 <- set_attributes(SVG2, 
      node = grep("marker", SVG2@summary$id, value = TRUE),
      node_attr = "id",
      attr = "transform",
      pattern = "scale\\(-0.2\\)",
      replacement = "scale(-0.15)")
    
    # export map as SVG again
    #write_svg(SVG2, file = paste0("../data/simulation/central_metabolism_", df$carbon_source[1], ".svg"))
  }
) %>% invisible()

Metabolic flux with CO2 refixation	Metabolic flux on fructose

Metabolic flux on succinate	Metabolic flux on formic acid