Constraints for an R. eutropha resource allocation model

1 Description

This R notebook is a bioinformatics pipeline to collect constraints for a genome scale, resource allocation model in the chemolithoautotroph Ralstonia eutropha (a.k.a. Cupriavidus necator).

A resource allocation model can be coarse-grained (few symbolic reactions) or have genome scale detail (all known biochemical reactions and their associated genes). However, both types of models need to be constrained by a set of parameters to make realistic predictions. Depending on the model frame work, constraints can be equality constraints (example: turnover number of an enzyme E kcat_E = 100 s^-1), or inequality constraints (0 s^-1 <= kcat_E <= 100 s^-1).

This notebook has the purpose to collect constant and growth-rate dependent constraints as they are used in RBA models. In RBApy, apparent enzyme efficiencies (k_app), protein abundance, molecular machine abundance (protein/macromolecule complexes), and fluxes can be constrained. RBApy has the following possibilities for custom constraints.

• constants (example: A = 0.1)
• linear relationship, e.g. with growth rate µ (example: B = 2 * µ + 0.1)
• Michaelis-Menthen like kinetics for k_app (example: kapp = kcat * [S] / ([S] + Km))

Different types of data were used to constrain the resource allocation model. The primary data is protein abundance determined by mass spectrometry for R. eutropha using different growth rates and carbon sources. This data is used to estimate and constrain k_app, enzyme abundance, and non-enzyme protein abundance.

2 Libraries

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

3 k_app estimation using proteomics/fluxomics

3.1 Background

One approach to estimate k_app values is to use enzyme abundance data and estimated or measured fluxes. The k_app is a parameter that links enzyme fluxes to enzyme abundance, so that flux v = k_app * [E]. The k_app includes also the saturation s of the enzyme, so that k_app ~ k_cat * s. The saturation is between 0 and 1 so that k_app is lower or equal to k_cat.

However, the idea is that the k_app parameter is kept constant for different conditions, while the flux (and the enzyme abundance) is allowed to change. We can then find different scenarios for enzyme saturation states:

• the flux through an enzyme and therefore its abundance increases in the model, but not in measured enzyme abundance: –> k_app is higher than estimated, enzyme not saturated
• the flux through an enzyme and therefore its abundance increases in the model, and in measured enzyme abundance: –> k_app is correctly estimated, enzyme is saturated

3.2 Fluxomics and proteomics input

To estimate k_app from available data, we only need to reformulate the simplified rate equation to: k_app = v / [E]. For this purpose, we need two types of information:

1. enzyme abundance in mmol gDCW^-1 for known substrate uptakes rates (chemostat cultivations, maximum growth rate 0.25 h^-1). The maximum tested growth rate was chosen to obtain realistic enzyme saturation. Protein abundance data was obtained by MS measurements. The relative protein mass fraction was determined by dividing MS intensity per protein by sum of all intensities. The mass fraction (g/g total protein) was then converted to protein concentration in mmol/gDCW by multiplying it with estimated protein concentration of 0.65 g protein/gDCW (done previously), and then dividing by molar mass of each protein (g/mol).
2. flux per enzyme, in mmol h^-1 gDCW^-1 for the corresponding substrate uptake rates under point 1. This information is obtained from flux sampling (FBA) simulations that were constrained to realistic flux distributions using data from Alagesan et al., 2017. Input and output fluxes (for exchange reactions) were determined from chemostat cultivations.

---

Determine enzyme abundance for standard conditions. Import Ralstonia mass fractions for the four tested conditions, fructose, succinate, formate, and ammonium (N-) limitation (with C-source fructose). Growth rate in all conditions was fixed to µ = 0.25 h^-1, the maximum growth rate used in chemostat experiments.

# import proteomics data
load("../data/input/Ralstonia_eutropha.Rdata")

# simplify condition strings
Ralstonia_eutropha <- Ralstonia_eutropha %>% ungroup %>%
  mutate(condition = sapply(condition, function(x){
    unlist(stri_extract_all_words(x))[1]}
  ))

# make new data frame to merge fluxomics and proteomics data
df_flux_protein <- Ralstonia_eutropha %>%
  
  # select only required columns and filter for highest mu
  mutate(condition = recode(condition, `FA` = "formate", `NLIM` = "ammonium", `FRC` = "fructose", `SUC` = "succinate")) %>%
  select(condition, locus_tag, growthrate, mass_g_per_gDCW, MolWeight) %>%
  filter(growthrate == 0.25) %>% 
  
  # calculate protein concentration in mmol/gDCW.
  # MW in kDa must be converted to g/mol, and concentration to mmol
  mutate(
    conc_mmol_gDCW = mass_g_per_gDCW * 1000 / (MolWeight * 1000),
    conc_mmol_gDCW = replace_na(conc_mmol_gDCW, 0)
  )

head(df_flux_protein)

# export to csv
for (cond in unique(df_flux_protein$condition)) {
  filter(df_flux_protein, condition == cond) %>% 
    select(-condition) %>%
    write_csv(paste0("../data/simulation/kapp_fitting/protein_concentration_", cond, ".csv"))
}

---

Determine reaction fluxes for same condition. Several approaches were tested to obtain minimal and maximal fluxes per enzyme. The best approach turned out to be flux sampling with the additional constraint of reaching at least 95% of the maximum growth rate found by FBA. Flux sampling with 100 iterations was performed using COBRApy. We can see that there is a threshold of standard deviation ~ 1 above which variation gets very high for some reactions. These are artificial cycles. The exception to this is N-limitation that can have higher flux variability due to the surmount carbon supply (fructose is not limiting).

# import flux sampling results for four conditions, max tested growth rate
df_sampling <- lapply(list.files("../data/simulation/kapp_fitting/", pattern = "^FSA", full.names = TRUE), 
  read_csv, col_types = cols()) %>% bind_rows(.id = "condition") %>%
  mutate(condition = recode(condition, `1` = "formate", `2` = "succinate", 
    `3` = "ammonium", `4` = "fructose")) %>%
  select(-2) %>% 
  gather(key = "reaction_id", value = "flux", -condition) %>%
  
  group_by(condition, reaction_id) %>% summarize(
    flux_mmol_gDCW_h = median(flux), 
    sd_flux = sd(flux), 
    min_flux = min(flux), 
    max_flux = max(flux),
    CV = abs(sd_flux/flux_mmol_gDCW_h)) %>%
  filter(!grepl("EX_", reaction_id)) %>%
  arrange(desc(sd_flux)) %>%
  
  # join with single optimal solution from loopless FBA.
  # Will help us to identify free running cycles in FSA
  left_join(
    lapply(list.files("../data/simulation/kapp_fitting/", pattern = "^FBA", full.names = TRUE),
    read_csv, col_types = cols()) %>% bind_rows(.id = "condition") %>%
    mutate(condition = recode(condition, `1` = "formate", `2` = "succinate", 
      `3` = "ammonium", `4` = "fructose")) %>%
    rename(reaction_id = X1, flux_mmol_gDCW_h_FBA = loopless)
  )

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

  
# plot results; 
plot_sampling <- lapply(list(A = c(-50, 550), B = c(-3, 28)), function(ylim) {
  xyplot(sd_flux ~ seq_along(sd_flux) | condition,
    filter(df_sampling, sd_flux != 0), 
    between = list(x = 0.5, y = 0.5),
    groups = condition, xlab = NULL, ylim = ylim,
    scales = list(alternating = FALSE),
    par.settings = custom.colorblind(), layout = c(4, 1),
    panel = function(x, y, ...) {
      panel.grid(h = -1, v = -1, col = grey(0.9))
      panel.abline(h = 1, lty = 2, col = grey(0.5))
      panel.xyplot(x, y, ...)
      panel.key(..., points = FALSE, corner = c(0.95,0.95))
    }
  )
})

print(plot_sampling[[1]], split = c(1,1,1,2), more = TRUE)
print(plot_sampling[[2]], split = c(1,2,1,2))

# process flux distribution, e.g. by removing extremely high fluxes
df_sampling2 <- df_sampling %>%
  
  # Filter out summary and outdated reactions
  filter(!reaction_id %in% c("Biomass", "Maintenance", "Protein", "Carbohydrate", 
    "Phospholipid", "   Cofactors_and_vitamins", "CBBCYC", "PYK1", "PYK2", "PYK3",
    "DHFR2", "DHFR3", "DHFR2p", "DHFR3p")) %>%
  
  # replace extreme fluxes with a min and max estimated from loopless FBA
  mutate(
    min_flux = if_else(sd_flux > 1, -abs(flux_mmol_gDCW_h_FBA), min_flux),
    max_flux = if_else(sd_flux > 1,  abs(flux_mmol_gDCW_h_FBA), max_flux)
  ) %>%
  
  # add an error margin to the sampled min and max fluxes, to help the solver
  # find a feasible solution
  mutate(
    min_flux = min_flux-(abs(min_flux)*1), 
    max_flux = max_flux+(abs(max_flux)*1)
  ) %>%
  
  # re-formatting of table
  arrange(desc(abs(max_flux))) %>%
  mutate(reaction_id = paste0("R_", reaction_id))

head(df_sampling2)

# export to csv
for (cond in unique(df_sampling2$condition)) {
  filter(df_sampling2, condition == cond) %>% ungroup %>%
    select(-condition) %>%
    write_csv(paste0("../data/simulation/kapp_fitting/model_flux_sampling_", cond, ".csv"))
}

3.3 Determine k_app values

The final step is to merge both datasets by computing the enzyme abundance allocated to each reaction. The estimated k_app is then determined by dividing the flux through the enzyme abundance available for this reaction. This step was not performed in this R notebook but using the RBA built-in functions from the RBApy estim folder. Briefly, flux boundaries and protein concentration data was exported and serves as input for the script kapp.py (link). This script performs a series of FBA and RBA simulations constrained by the input fluxes and protein concentrations. It then gives an estimation of the k_app for a particular condition.

4 Fraction of protein per compartment

The RBA model takes as another input parameter (or constraint) the fraction of protein per compartment. This constraint is important as it allows the cell to have only a limited amount of protein in cytoplasm or membrane compartments, for example. This constraint can be constant or it can be growth rate dependent e.g. by a linear relationship.

The first step is to prepare a table with protein abundance and localization. Protein abundance can be in any unit according to the RBApy manual, but it’s best to use mol fraction instead of mass fraction, as all other RBApy functions also use mmol. The mol fraction is already available in the processed data set. The built-in estim functions don’t seem to be well supported in python3 and raise errors. We therefore do the estimation manually by fitting linear models to the growth rate-protein abundance relationship. The idea is similar to the RBA estim function but less complicated. We focus on the standard condition as a training data set (fructose as carbon source, no NH4+ limitation).

df_protein <- Ralstonia_eutropha %>% 
  
  # select only required columns and spread to long format
  select(condition, locus_tag, growthrate, Psortb_localization, mol_fraction) %>%
  set_names(c("condition", "protein", "growthrate", "location", "mol_fraction")) %>%
  mutate(condition = recode(condition, `FA` = "formate", `NLIM` = "ammonium", `FRC` = "fructose", `SUC` = "succinate")) %>%
  
  # match localization names to model, simplify
  mutate(location = recode(location, Unknown = "Cytoplasm", Cytoplasmic = "Cytoplasm")) %>%
  mutate(location = replace(location, location != "Cytoplasm", "Cell_membrane"))

head(df_protein)

Now we can summarize the data by taking the sum of mol fraction over condition and localization. A simple approach to finding linear functions, where all proteins of all locations sum to one for a specific growth rate, would be to fit linear models for all compartments except one (e.g. cytoplasmic proteins, the biggest compartment). This one will then get a linear model fitted from the residual protein mass. We would expect the error for Cytoplasm to be negligibly small as it is the biggest compartment. However it turned out that the linear models fitted to both compartments perfectly sum to one already (see below). An estimation from residual protein fraction is therefore not necessary.

# First retrieve the proteins that are not part of the model, which is needed
# to  calculate non-enzymatic fraction later on
model_genes <- c(
  read_csv("../data/input/model_reactions.csv", col_types = cols()) %>% separate_rows(genes, sep = ", ") %>%
    filter(!duplicated(genes)) %>% pull(genes),
  read_tsv("../data/simulation/macro_machines/ribosome.tsv", col_types = cols())[["Gene names"]],
  read_tsv("../data/simulation/macro_machines/chaperones.tsv", col_types = cols())[["Gene names"]],
  read_tsv("../data/simulation/macro_machines/transcription.tsv", col_types = cols())[["Gene names"]],
  read_tsv("../data/simulation/macro_machines/replication.tsv", col_types = cols())[["Gene names"]]
)

Warning: Missing column names filled in: 'X1' [1]

# extract locus_tags only
model_genes <- model_genes %>% stri_extract_first_regex(
  pattern = "H16_[A-Z][0-9]{4}|PHG[0-9]{3}")

# add new NE_protein column
df_prot_per_comp <- df_protein %>%
  mutate(NE_protein = !protein %in% model_genes) %>%
  group_by(condition, location, growthrate)

The following plot shows a slight increase in cytoplasmic proteins and decrease in cell membrane proteins with growth rate.

plot_prot_comp <- xyplot(prot_per_compartment ~ growthrate, 
  df_prot_per_comp %>%
    summarize(prot_per_compartment = sum(mol_fraction, na.rm = TRUE)),
  groups = location, ylim = c(-0.1, 1.1),
  par.settings = custom.colorblind(), cex = 0.7,
  scales = list(alternating = FALSE),
  between = list(x = 0.5, y = 0.5),
  xlab = expression("µ [h"^-1*"]"),
  ylab = "protein mass fraction",
  panel = function(x, y, ...) {
    panel.grid(h = -1, v = -1, col = grey(0.9))
    panel.superpose(x, y, ...)
    panel.key(corner = c(0.1, 0.55), ...)
  }, panel.groups = function(x, y, ...) {
    panel.xyplot(x, y, ...)
    panel.lmlineq(x, y, r.squared = TRUE,
      pos = 3, offset = 1, ...)
  }
)

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

print(plot_prot_comp)

5 Fraction of non-enzymatic protein per compartment

The previous calculation determined how the total protein pool is distributed over compartments. In this section, we will go one level deeper and determine the fraction of non-enzymatic protein per compartment, analogously to the previous calculation. That means every compartment’s protein pool is further subdivided into two sectors, enzymatic and non-enzymatic proteins, and these sectors can have –again– a linear dependency on growth rate. According to the RBApy manual, the proteins that are considered non-enzymatic are all proteins not acting as enzymes or molecular machines in the model (such as ribosomes, chaperones, DNA polymerase and so on). Non-enzymatic (NE) proteins are therefore all proteins not captured by the RBA model at all. The fraction of NE proteins per compartment is a percentage of the proteins for that compartment only (each compartment sums to 1).

# determine NE protein per compartment
df_ne_prot_per_comp <- df_prot_per_comp %>% group_by(NE_protein, .add = TRUE) %>%
  summarize(prot_per_compartment = sum(mol_fraction, na.rm = TRUE)) %>%
  summarize(ne_prot_per_compartment = 
    prot_per_compartment[2]/sum(prot_per_compartment)
  )

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

  
# the fraction of NE protein per compartment decreases with growth rate
plot_prot_NE <- xyplot(ne_prot_per_compartment ~ growthrate,
  df_ne_prot_per_comp,
  groups = location, ylim = c(-0.1, 1.1),
  par.settings = custom.colorblind(),
  scales = list(alternating = FALSE),
  between = list(x = 0.5, y = 0.5), cex = 0.7,
  xlab = expression("µ [h"^-1*"]"),
  ylab = "NE protein mass fraction",
  panel = function(x, y, ...) {
    panel.grid(h = -1, v = -1, col = grey(0.9))
    panel.superpose(x, y, ...)
    panel.key(corner = c(0.1, 0.1), ...)
  }, panel.groups = function(x, y, ...) {
    panel.xyplot(x, y, ...)
    panel.lmlineq(x, y, r.squared = TRUE,
      pos = 3, offset = 1, ...)
  }
)

print(plot_prot_NE)

---

Summary: Regardless of the cultivation condition or substrate limitation, it is possible to see growth rate dependent trends. Overall, the fraction of cytoplasmic proteins increases slightly with growth rate while membrane-associated proteins decrease. This makes sense as often expression of transporters is down-regulated if bacteria don’t need to scavenge substrates. We also see that the absolute majority of proteins is located in the cytoplasm (> 85% under all conditions).

The non-enzymatic fraction of proteins (proteins not covered by the model) decreases with growth rate from around 60% o 50% for cytoplasmic proteins, but is relatively constant for membrane associated ones. Most membrane-associated proteins (90% of mol fraction) are not covered by the model. The following figure summarizes all growth rate dependent trends.

plot_prot_pie <- df_prot_per_comp %>% group_by(NE_protein, .add = TRUE) %>%
  summarize(prot_per_compartment = sum(mol_fraction, na.rm = TRUE)) %>%
  ungroup %>% mutate(
    location = recode(location, Cytoplasm = "CP", Cell_membrane = "CM"),
    location = paste(location, ifelse(NE_protein, "NE", "E"), sep ="-")) %>%
  filter(growthrate == 0.25) %>%
  #filter(condition == "FRC") %>%
  
  xyplot( ~ prot_per_compartment | condition, .,
    groups = factor(location), scales = list(draw = FALSE),
    xlab = "mol fraction of proteome, per compartment",
    par.settings = custom.colorblind(), layout = c(4,1),
    between = list(x = 0.5, y = 0.5),
    cex = 0.6, border = grey(0.3),
    panel = function(x, y, ...) {
      panel.piechart(x, diameter_inner = 0.1, 
        diameter_sector = 0.15, ...)
    }
  )

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

print(plot_prot_pie)

---

6 Combining estimated k_app values to global consensus k_app

The purpose of this section is to merge different k_app estimations from different conditions to obtain a consensus k_app estimation. The reasoning here is that k_app is estimated from dividing the (sampled or measured) metabolic flux through an enzyme by its abundance, providing an estimate of enzyme efficiency. This efficiency depends also on the saturation of the enzyme. Estimation of k_app was therefore performed for highest available growth rates/fluxes for four different conditions, in order to obtain enzyme efficiency at the highest saturation among the tested conditions.

Example: Formate dehydrogenase is not used under growth on fructose (no or low efficiency) while it is used under growth on formate (high efficiency). The maximum k_app should therefore be taken from formate condition.

6.1 Import k_app estimation data obtained from RBA estim

# data directory
kapp_calib_dir <- "../data/simulation/kapp_fitting/"

# setting global parameter
max_saturation = 0.75

# load kapp estimation from different conditions
conditions <- c("ammonium", "formate", "fructose", "succinate")
kapp_files <- paste0(kapp_calib_dir, "kapp_estimate_", conditions, ".csv")

df_kapp <- lapply(kapp_files, function(df) {
  read_tsv(df, col_names = c("reaction_id", "kapp_forward", "kapp_reverse"))}) %>% 
  bind_rows(.id = "condition") %>%
  
  # change condition names
  mutate(condition = recode(condition, !!!(conditions %>% setNames(1:4)))) %>%
  
  # add a saturation factor S, increasing kapp by 1/S
  mutate_at(vars(contains("kapp_")), function(x) x/max_saturation) %>%
  
  # apply median normalization to account for different in silico growth 
  # rates in FBA simulation (different total flux)
  group_by(condition) %>% 
  mutate(kapp_forward = kapp_forward*(median(.[["kapp_forward"]])/median(kapp_forward))) %>%
  mutate(kapp_reverse = kapp_forward) %>%
  
  # sort decreasingly by average kapp
  group_by(reaction_id) %>% 
  mutate(kapp_mean = mean(kapp_forward)) %>%
  arrange(desc(kapp_mean))

── Column specification ─────────────────────────────────────────────────────────────────────────────────────────────
cols(
  reaction_id = col_character(),
  kapp_forward = col_double(),
  kapp_reverse = col_double()
)


── Column specification ─────────────────────────────────────────────────────────────────────────────────────────────
cols(
  reaction_id = col_character(),
  kapp_forward = col_double(),
  kapp_reverse = col_double()
)


── Column specification ─────────────────────────────────────────────────────────────────────────────────────────────
cols(
  reaction_id = col_character(),
  kapp_forward = col_double(),
  kapp_reverse = col_double()
)


── Column specification ─────────────────────────────────────────────────────────────────────────────────────────────
cols(
  reaction_id = col_character(),
  kapp_forward = col_double(),
  kapp_reverse = col_double()
)

head(df_kapp)

# plot kapp value distribution per condition
plot_kapp_dist <- lapply(c("condition", "all_conditions"), function(cond) {
  xyplot(log10(kapp_forward) ~ as.numeric(factor(reaction_id, unique(reaction_id))) | get(cond), 
    df_kapp %>% mutate(all_conditions = "all conditions"), 
    groups = condition, cex = 0.7, pch = 19,
    xlab = "number of reactions", ylab = expression("log"[10]*"  k"[app]*" [h"^-1*"]"),
    as.table = TRUE, between = list(x = 0.5, y = 0.5),
    par.settings = custom.colorblind(),
    scales =list(alternating = FALSE),
    panel = function(x, y, ...){
      panel.grid(h = -1, v = -1, col = grey(0.9))
      panel.xyplot(x, y, ...)
      med = median(y)
      panel.abline(h = med, lty = 2, lwd = 2, col = grey(0.5))
      if (cond == "all_conditions") {
        panel.text(grid::unit(175, "npc"), med+3,
          labels = paste0("median kapp = ", round(10^med), " [h-1]"))
        panel.key(..., corner = c(0.98, 0.95))
      }
    }
  )
})

print(plot_kapp_dist[[2]])

6.2 Manual adjustment of k_app for selected reactions

We import a table with correction factors for k_app values that have been estimated incorrectly. This set of k_app values was determined by comparing the correlation between predicted and measured protein abundance. Proteins that are more than 10 fold away from measured abundance are candidates for manual correction. Other reactions that need correction are missing k_app values for important transporters. These are also included in the following section.

# set of adjusted kapp values for outlier enzymes not covered in kapp estimation
kapp_median <- median(df_kapp$kapp_forward)
kapp_extra <- read_csv(paste0(kapp_calib_dir, "kapp_extra.csv"))

── Column specification ─────────────────────────────────────────────────────────────────────────────────────────────
cols(
  reaction_id = col_character(),
  ratio_mol_fraction = col_double()
)

# print mean and sd of log-transformed kapp distribution
df_kapp$kapp_forward %>% log10 %>% mean

[1] 3.813431

df_kapp$kapp_forward %>% log10 %>% sd

[1] 1.115881

# kapp estimation for reactions where information is missing
kapp_missing <- kapp_extra %>%
  filter(!(reaction_id %in% df_kapp$reaction_id)) %>%
  mutate(
    kapp_forward = ratio_mol_fraction*kapp_median,
    kapp_reverse = kapp_forward) %>%
  select(-ratio_mol_fraction)

# kapp estimation for transporters
kapp_transport <- data.frame(stringsAsFactors = FALSE,
  reaction_id = c("R_H2Ot_enzyme", "R_CO2t_enzyme", "R_O2t_emzyme", "R_SO4t_enzyme",
    "R_PIt2r_enzyme", "R_Ktr_enzyme", "R_NAt3_15_enzyme", "R_FE2abc_enzyme",
    "R_MG2t_enzyme", "R_COBALTt5_enzyme", "R_FORt_enzyme", "R_NH4t_enzyme", 
    "R_FRUabc_enzyme", "R_SUCCt2r_enzyme", "R_FRUpts2_enzyme", 
    "R_SUCCabc_enzyme", "R_SUCFUMt_enzyme"),
  kapp_forward = c(rep(1000000, 11), rep(100000, 3), rep(0, 3))
  ) %>% mutate(kapp_reverse = kapp_forward)

The final consensus kapp table is prepared by adding missing kapps for transport reactions and correcting several falsely estimated kapps. The table is exported as *.csv file. RBApy requires the table in a specific format:

• tab-separated values, no header
• units in 1/h instead of 1/s (x 3600)
• reaction ID in column 1, format (R_)ID_enzyme or ID_transporter
• max k_kat in column 2
• min k_kat in column 3 (backward efficiency, can be identical)
• no NA values (only complete rows)

# Select maximum as aggregation metric
df_kapp_export <- df_kapp %>% group_by(reaction_id) %>%
  summarize(
    kapp_forward = max(kapp_forward),
    kapp_reverse = max(kapp_reverse)
  ) %>%
  
  # add missing kapps
  filter(!(reaction_id %in% c("R_HEXf_enzyme","R_FRUabc_enzyme"))) %>%
  bind_rows(kapp_missing) %>%
  bind_rows(kapp_transport) %>%
  
  # update kapp estimation for some reactions where estimation was inaccurate
  left_join(filter(kapp_extra, reaction_id %in% df_kapp$reaction_id)) %>%
  mutate(ratio_mol_fraction = replace_na(ratio_mol_fraction, 1)) %>%
  mutate(
    kapp_forward = kapp_forward * ratio_mol_fraction,
    kapp_reverse = kapp_reverse * ratio_mol_fraction
  ) %>%
  select(-ratio_mol_fraction)

Joining, by = "reaction_id"

# export final result
write_tsv(df_kapp_export, "../data/simulation/kapp_fitting/kapp_consensus.csv", col_names = FALSE)

Finally we can create a summary figure for the parameters constraining the RBA model. The first is k_app estimation and the two other constraints are proteome mass fraction per compartment and proteome mass fraction of enzymatic (all RBA model reactions) and non-enzymatic proteins, per compartment.

print(plot_kapp_dist[[2]], position = c(0, 0.51, 1, 1), more = TRUE)
print(plot_prot_comp, position = c(0, 0, 0.5, 0.55), more = TRUE)

print(plot_prot_NE, position = c(0.5, 0, 1, 0.55), more = TRUE)
grid::grid.text(c("A", "B", "C"), x = c(0.02, 0.02, 0.5), 
  y = c(0.97, 0.5, 0.5))