Physiological measurements from R. eutropha chemostats

1 Description

This R notebook contains data and analysis regarding physiological measurements obtained from chemostat cultivations of Ralstonia eutropha (a.k.a. Cupriavidus necator).

2 Libraries

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

3 Chemostat cultivation, biomass and growth rate

3.1 Import and reshape cultivation data

First we load the three cultivation data tables. These contain all four limiting conditions. We are thinning out the OD measurement data by 75% so that file sizes become smaller and more plottable.

# load chemostat measurements
df_chem <- read.csv("../data/input/20180921_chemostat_OD.csv") %>%
  
  # remove outlier data points originating from sampling
  filter(!od_value > 0.6) %>%
  
  # thinning data to reduce file size
  slice_sample(n = 4000)

3.2 OD and growth rate

The next step is to plot the OD₇₂₀ and growth rate for all cultivations.

plot_OD <- xyplot(od_value ~ batchtime_h | condition, 
  arrange(df_chem, batchtime_h),
  groups = replicate, par.settings = custom.colorblind(),
  layout = c(1,4), as.table = TRUE,
  type = "l", xlim = c(0, 335),
  cols = custom.colorblind()$superpose.polygon$col[1:4],
  between = list(x = 0.5, y = 0.5), lwd = 2,
  ylim = c(-0.06, 0.86), ylab = expression("OD"[720*nm]),
  panel = function(x, y, cols, ...) {
    panel.xyplot(x, y, col = cols[panel.number()], ...)
  }
)

plot_dil <- xyplot(dilution_rate ~ batchtime_h | condition, 
  filter(df_chem, replicate == 1), 
  par.settings = custom.colorblind(), type = "l",
  between = list(x = 0.5, y = 0.5), ylim = c(-0.03, 0.43),
  xlab = "time [h]", ylab = expression("µ [h"^-1*"]"),
  panel = function(x, y, ...){
    panel.grid(h = -1, v = -1, col = grey(0.9))
    panel.horizonplot(c(x, 335), c(y, max(y)), col.regions = grey(0.7, alpha = 0.5))
    arrow_x <- c(sort(x)[which(diff(sort(y)) != 0)], max(x)-1)
    panel.arrows(x0 = arrow_x, y0 = 0.4, x1 = arrow_x, y1 = 0.32,
      type = "closed", col = grey(0.5), lwd = 2, fill = grey(0.5), 
      length = 0.05, unit = "npc")
  }
)

print(doubleYScale(plot_OD, plot_dil, use.style = FALSE, add.ylab2 = TRUE, under = TRUE))

3.3 Residual substrate concentration

The residual substrate (carbon sources only) was determined using HPLC measurement. In order to merge cultivation and residual substrate measurements, we summarize average OD and growth rate for all cultivations leaving out the first 25 h batch growth phase.

df_summary <- filter(df_chem, batchtime_h > 25) %>%
  
  # summarize average OD
  group_by(condition, dilution_rate, replicate) %>%
  summarize(OD720 = mean(od_value, na.rm = TRUE)) %>%
  
  # then we need medium substrate concentration for FA, FRC, NH4CL, SUC
  mutate(c_substrate = case_when(
    condition == "ammonium" ~ 0.025,
    condition == "formate" ~ 1.0,
    condition == "fructose" ~ 0.5,
    condition == "succinate" ~ 0.5
  ))

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

The residual substrate measurement is imported, and merged with the cultivation summary table. For concnetration of substrate that is really taken up is calculated by subtracting the residual substrate (mg) from feed substrate concentration (g).

df_residual <- read_csv("../data/input/20181116_chemostat_residual_substrate.csv", col_types = cols()) %>%
  rename(condition = limitation) %>%
  mutate(condition = recode(condition, `formic acid` = "formate", `Nlim - fructose` = "ammonium")) %>%
  mutate(c_residual_g_L = mean_concentration_mg.L/1000 %>% replace(., . < 0, 0))

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

# join the summary df with residual substrate data
df_summary <- left_join(df_summary, df_residual %>%
    select(condition, dilution_rate, c_residual_g_L) %>%
    filter(!duplicated(c_residual_g_L), condition != "ammonium")) %>%
  mutate(c_residual_g_L = c_residual_g_L %>% replace_na(0))

Joining, by = c("condition", "dilution_rate")

print(head(df_summary))

plot_residual <- xyplot(concentration_mg.L ~ factor(dilution_rate) | condition,
  df_residual,
  par.settings = custom.colorblind(), as.table = TRUE,
  pch = 19, ylim = c(-50, 1450), layout = c(4, 1),
  between = list(x = 0.5, y = 0.5),
  scales = list(alternating = FALSE, x = list(at = c(1,3,5))),
  xlab = expression("µ [h"^-1*"]"), ylab = "concentration [mg/L]",
  panel = function(x, y, ...) {
    panel.grid(h = -1, v = -1, col = grey(0.9))
    panel.barplot(x, y, ewidth = 0.3, lwd = 2, fill = "white", fill_alpha = 1, ...)
    panel.xyplot(x, y, alpha = 0.3, ...)
  }
)

print(plot_residual)

3.4 DCW and PHB content of cells

The next step is to add biomass concentration (DCW) in g/L. For carbon limitations, the DCW was determined in a separate series of 50 mL shake flask cultivations with substrate concentration identical to chemostat experiments (fructose: 0.5 g/L, succinate 0.5 g/L, 4 replicates). The concentration of formate was doubled to obtain sufficient biomass for DCW measurement (2.0 g/L). For ammonium limitation, the cultivation was reproduced in chemostats with an N-limiting concentration of 0.05 g/L NH4Cl and non-limiting carbon concentration of 2.0 g/L fructose.

df_DCW_PHB <- bind_rows(
  read_csv("../data/input/20210810_gDCW_carbon_lim.csv", col_types = cols()),
  read_csv("../data/input/20210810_gDCW_nitrogen_lim.csv", col_types = cols())
) %>%
  filter(condition != "ammonium WT") %>%
  mutate(condition = condition %>% str_replace(" ", "\n") %>%
    str_replace("ammonium", "ammon."))

# plot total gDCW and total gPHB
plot_PHB_Clim <- xyplot(g_DCW_L + g_PHB_L ~ factor(condition),
  filter(df_DCW_PHB, !str_detect(condition, "ammon")),
  par.settings = custom.colorblind(), ylim = c(0, 0.3),
  lwd = 2, xlab = "", ylab = expression("g L"^-1),
  panel=function(x, y, ...) {
    panel.grid(h = -1, v = -1, col = grey(0.9))
    panel.barplot(x, y, beside = TRUE, draw_points = TRUE, ...)
    panel.annotate(x, y, offset = 0.04, digits = 3, beside = TRUE, cex = 0.7, ...)
    panel.key(groups = c("g DCW/L", "g PHB/L"), corner = c(0.95, 0.95))
  }
)

plot_PHB_Nlim <- xyplot(g_DCW_L + g_PHB_L ~ factor(condition, unique(condition)[c(2:6,1)]),
  filter(df_DCW_PHB, str_detect(condition, "ammon")),
  par.settings = custom.colorblind(), ylim = c(0, 0.3),
  lwd = 2, xlab = "", ylab = expression("g L"^-1),
  panel=function(x, y, ...) {
    panel.grid(h = -1, v = -1, col = grey(0.9))
    panel.barplot(x, y, beside = TRUE, draw_points = TRUE, ...)
    panel.annotate(x, y, offset = 0.04, digits = 3, beside = TRUE, cex = 0.7, ...)
    panel.key(groups = c("g DCW/L", "g PHB/L"), corner = c(0.95, 0.95))
  }
)

# plot yield on substrate (excluding PHB)
plot_yield_Clim <- xyplot(g_DCW_g_substrate_noPHB ~ factor(condition),
  filter(df_DCW_PHB, !str_detect(condition, "ammon")),
  par.settings = custom.colorblind(), ylim = c(-0.02, 0.52),
  lwd = 2, xlab = "", ylab = expression("g DCW g S"^-1),
  panel=function(x, y, ...) {
    panel.grid(h = -1, v = -1, col = grey(0.9))
    panel.barplot(x, y, draw_points = TRUE, ...)
    panel.annotate(x, y, offset = 0.1, digits = 3, cex = 0.7, ...)
  }
)

plot_yield_Nlim <- xyplot(g_DCW_g_substrate_noPHB ~ 
    factor(condition, unique(condition)[c(2:6,1)]),
  filter(df_DCW_PHB, str_detect(condition, "ammon")),
  par.settings = custom.colorblind(), ylim = c(-0.1, 2.4),
  lwd = 2, xlab = "", ylab = expression("g DCW g S"^-1),
  panel=function(x, y, ...) {
    panel.grid(h = -1, v = -1, col = grey(0.9))
    panel.barplot(x, y, draw_points = TRUE, ...)
    panel.annotate(x, y, offset = 0.5, digits = 3, cex = 0.7, ...)
  }
)

print(plot_PHB_Clim, position = c(0,0.49,0.42,1), more = TRUE)
print(plot_PHB_Nlim, position = c(0.38,0.463,1,1), more = TRUE)

print(plot_yield_Clim, position = c(0.01,0,0.42,0.51), more = TRUE)
print(plot_yield_Nlim, position = c(0.39,-0.027,1,0.51))

grid::grid.text(c("A", "B", "C", "D"), x = c(0.03, 0.4, 0.03, 0.4),
  y = c(0.97, 0.97, 0.48, 0.48), gp = grid::gpar(cex = 1.2))

3.5 Calculate biomass yield and maintenance

The yield on formate was considerably lower than the maximum biomass yield on formate (0.093 g DCW / g S, Grunwald et al., 2015). The optical density of the cultures was 40% lower than the OD in bioreactors. The explanation is the inhibitory effect of formate. The yield was determined with a starting concentration of 2 g/L which was shown to reduce yield significantly. Grunwald et al. established a linear model to estimate the toxic effect of residual formic (RF) acid concentration: fraction of max yield = (0.169 - 0.049 * RF)/0.169. For RF = 2.0 g/L, the obtained biomass yield is only 0.42 of theoretical maximum. In order to obtain a realistic estimation of yield in chemostats with no or low residual formate, formate yield was corrected based on above relationship for the purpose of modeling.

The substrate uptake rate qS is calculated using the corrected substrate concentration (S_uptake = S_feed - S_residual). The rate qS is then dilution rate times feed concentration divided by experimentally determined biomass in the culture (excluding PHB): qS [g/gDCW] = D [h-1] * c [g/L] / biomass [gDCW/L]. The DCW for formate and ammonium was determined using twice the original concentration, hence the g DCW is divided by two for compatibility.

df_summary <- df_DCW_PHB %>%
  filter(condition != "ammon.\nPHB-") %>%
  mutate(condition = str_replace(condition, "ammon.\nu=.*", "ammonium")) %>%
  group_by(condition, growth_rate) %>%
  
  # summarize replicates to mean biomass and PHB
  summarise(.groups = "drop",
    g_DCW_L = mean(g_DCW_L_noPHB),
    g_PHB_L = mean(g_PHB_L)) %>%
  mutate(
    g_DCW_L = case_when(
      condition == "formate" ~ g_DCW_L/(0.42*2),
      condition == "ammonium" ~ g_DCW_L/2, TRUE ~ g_DCW_L),
    g_PHB_L = case_when(condition %in% c("formate", "ammonium") ~ g_PHB_L/2, TRUE ~ g_PHB_L),
    percent_PHB = g_PHB_L/(g_PHB_L+g_DCW_L)*100
  ) %>%
  
  # fill missing conditions (growth rates), and merge with main df
  rename(dilution_rate = growth_rate) %>%
  replace_na(list(dilution_rate = 0.05)) %>%
  complete(condition, dilution_rate) %>%
  fill(g_DCW_L, g_PHB_L, percent_PHB) %>%
  right_join(df_summary) %>%

  # substrate uptake rate per gDCW
  # qS [g/gDCW] = D [h-1] * c [g/L] / gDCW * L
  mutate(qS_g_gDCW = dilution_rate*(c_substrate-c_residual_g_L)/g_DCW_L) %>%
  
  # substrate uptake rate in C-mol per C-mol biomass (25.35 gDCW/Cmol = 0.03944 Cmol/gDCW)
  # (source: Grunwald et al., Mirco Biotech, 2015)
  mutate(Cmol_g = case_when(
    condition == "formate" ~ 1/46.03 * 1,
    condition == "fructose" ~ 1/180.16 * 6,
    condition == "succinate" ~ 1/118.09 * 4,
    condition == "ammonium" ~ 1/53.49 * 1 # here it is N-mol! NH4Cl has 1 N
  )) %>%
  
  # for NH4Cl, we actually use N-mol per N-mol biomass per h. 
  # N-mol in biomass is around 4x lower than C mol (bionumbers)
  mutate(qS_Cmol_Cmol_DCW = case_when(
    condition == "ammonium" ~ qS_g_gDCW * Cmol_g / (0.03944/4),
    TRUE ~ qS_g_gDCW * Cmol_g / 0.03944
  ))

Joining, by = c("condition", "dilution_rate")

Finally we can create a so-called Herbert-Pirt plot, that is growth rate versus substrate uptake rate. This plot will usually reveal a linear relationship, meaning that substrate uptake and growth are proportional to each other. If there is a kink (change in slope) of the data and one linear model would not be sufficient to describe the data, then the yield has changed, i.e. cells start using a different pathway for energy generation (e.g. respiration versus fermentation). In the default scenario (linear relationship), the yield is the slope of the linear model and the maintenance coeffcient is the intersection with the Y axis.

HP_plot <- function(data, xlimits = c(0, 0.3), ylimits = c(0, 1)) {
  xyplot(qS_g_gDCW ~ dilution_rate | condition, data,
    par.settings = custom.colorblind(), lwd = 1.5,
    xlim = xlimits, ylim = ylimits,
    ylab = expression("q"[S]*" [g h"^-1*" gDCW"^-1*"]"),
    xlab = expression('µ [h'^'-1'*']'),
    panel = function(x, y, ...) {
      panel.grid(h = -1, v = -1, col = grey(0.9))
      panel.xyplot(x, y, cex = 0.9, ...)
      # regression line through linear part of the data (omit µ = 0.25/h)
      if (data$condition[1] == "formate") sel <- x < 0.2 else sel <- x < 0.25
      x = x[sel]; y = y[sel]
      panel.lmline(x, y, ...)
      coef <- lm(y ~ x, data.frame(x, y))$coeff
      # displaying maintenance and yield coefficients
      panel.abline(h = coef[[1]], lty = 2, ...)
      panel.text(0.15, coef[[1]],
        paste("ms =", round(coef[[1]], 3), "g h-1 gDCW-1"),
        col = grey(0.3), pos = 3, cex = 0.7)
      panel.text(0.15, coef[[1]], paste(expression("Yx/S ="),
        round(1/coef[[2]], 3), "gDCW g-1"),
        col = grey(0.3), pos = 1, cex = 0.7)
    }
  )
}

print(HP_plot(data = filter(df_summary, condition == "formate"),
  ylimits = c(-7/4, 7)), position = c(0.017,0.5,0.5,1), more = TRUE)
print(HP_plot(data = filter(df_summary, condition == "fructose"),
  ylimits = c(-1.5/4, 1.5)), position = c(0.5,0.5,1,1), more = TRUE)

print(HP_plot(data = filter(df_summary, condition == "ammonium"),
  ylimits = c(-0.25/4, 0.25)), position = c(-0.017,0,0.5,0.5), more = TRUE)
print(HP_plot(data = filter(df_summary, condition == "succinate"),
  ylimits = c(-1.5/4, 1.5)), position = c(0.5,0,1,0.5))

---

For metabolic modeling, we shall also calculate the specific productivity of PHB as rate in g PHB per g DCW per hour. We can convert the unit to the convention used in metabolic modeling, mmol per g DCW per hour, using the molecular mass of the 3-hydroxy butyrate monomer, M = 104.1 g/mol, subtracted by one molecule water lost during the condensation into the polymer (M = 104.1 - 18 = 86.1 g/mol = 0.0861 g/mmol).

df_summary <- df_summary %>%
  mutate(
    mmol_PHB_g_DCW = percent_PHB/(100*0.0861),
    P_g_PHB_g_DCW_h = g_PHB_L/g_DCW_L*dilution_rate,
    P_mmol_PHB_g_DCW_h = P_g_PHB_g_DCW_h/0.0861
  )
  
xyplot(P_mmol_PHB_g_DCW_h ~ factor(dilution_rate) | condition, df_summary,
  par.settings = custom.colorblind(), layout = c(4, 1),
  as.table = TRUE, between = list(x = 0.5, y = 0.5),
  ylab = expression("P"[PHB]*" [mmol h"^-1*" gDCW"^-1*"]"),
  xlab = expression('µ [h'^'-1'*']'),
  scales = list(alternating = FALSE, x = list(at = c(1, 3, 5))),
  panel = function(x, y, ...) {
    panel.grid(h = -1, v = -1, col = grey(0.9))
    panel.xyplot(x, y, cex = 0.9, ...)
  }
)

# preview table with mass % and flux to PHB for all conditions
df_summary %>%
  select(condition, dilution_rate, percent_PHB, mmol_PHB_g_DCW, 
    P_g_PHB_g_DCW_h, P_mmol_PHB_g_DCW_h) %>%
  distinct %>% print

3.6 O₂ consumption and CO₂ emission

Outlet gas composition was analyzed with specific optical sensors for O₂ and CO₂. The gas uptake rate can then be determined knowing the input concentration, output concentration, flow rate per reactor, and biomass concentration per reactor. First load biomass and gas concentration raw measurements.

df_chemostat_Nlim <- read_csv("../data/input/20190716_chemostat_PHB.csv", col_types = cols()) %>%
  filter(mu <= 0.3) %>% mutate(condition = recode(
    condition, "C-limitation" = "fructose", "N-limitation" = "ammonium")) %>%
  mutate(replicate = if_else(channel > 4, channel-4, channel))

df_gas <- read_csv("../data/input/20190728_chemostat_gas.csv", col_types = cols())

Next we summarize all measurements to one measurement per condition and replicate. Then we use measurements for the replicates to calculate the mean gas uptake and emission rate normalized to biomass concentration.

df_gas <- df_gas %>%
  
  # remove mu higher than 0.3
  filter(mu <= 0.3) %>%
  
  # determine average delta air flow per cond and replicate
  group_by(condition, strain, mu, type, channel, `flow_rate [ml/min]`) %>%
  summarize(flow_mL_min = mean(`flow [ml/min]`)) %>%
  
  # renumber replicates from 1 to n
  group_by(condition, strain, mu, type) %>%
  mutate(replicate = seq_along(channel)) %>%
  arrange(condition, strain, mu, type, replicate) %>%
  
  # obtain biomass concentration from first table
  ungroup %>% left_join(
    group_by(df_chemostat_Nlim, strain, condition, mu) %>%
    summarize(gDCW_mL = mean(`DCW [g/ml]`, na.rm = TRUE)) %>% 
    ungroup
    )

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

Now that we have the biomass concentration in gDCW/mL, we can determine the CO₂ emission and O₂. The culture volume is 65 mL and the mass per volume CO₂ is 1.842 mg/mL (air at 1 atm and 20 °C). The mass per volume O₂ is 1.314 mg/mL (air at 1 atm and 20 °C) The CO₂ emission rate can then be calculated as:

qCO2 [mL min^-1 gDCW^-1] = flow_CO2 [mL min^-1] / (biomass gDCW mL^-1 * V_culture mL) qCO2 [g h^-1 gDCW^-1] =

df_gas <- df_gas %>% 
  
  # add volume to mass conversion
  mutate(m_gas_mg_mL = if_else(type == "CO2", 1.842, 1.314)) %>%
  mutate(
    q_mL_min_gDCW = flow_mL_min / (gDCW_mL * 65),
    q_g_h_gDCW = q_mL_min_gDCW * m_gas_mg_mL * 60 / 1000
  )

head(df_gas)

plot_CO2_em <- xyplot(q_g_h_gDCW ~ factor(mu) | condition,
  filter(df_gas, type == "CO2"), groups = strain,
  par.settings = custom.colorblind(), 
  pch = 19, lwd = 2, ylim = c(0, 1.5),
  as.table = TRUE, between = list(x = 0.5, y = 0.5),
  xlab = expression("µ [h"^-1*"]"), 
  ylab = expression("q"[CO2]*" g gDCW"^-1*" h"^-1),
  scales = list(alternating = FALSE),
  panel=function(x, y, ...) {
    panel.grid(h = -1, v = -1, col = grey(0.9))
    panel.errbars(x, y, ewidth = 0, ...)
    panel.key(..., which.panel = 2, corner = c(0.9, 0.9))
  }
)

plot_O2_up <- xyplot(q_g_h_gDCW*-1 ~ factor(mu) | condition,
  filter(df_gas, type == "O2"), groups = strain,
  par.settings = custom.colorblind(), 
  pch = 19, lwd = 2, ylim = c(0, 1.5),
  as.table = TRUE, between = list(x = 0.5, y = 0.5),
  xlab = expression("µ [h"^-1*"]"), 
  ylab = expression("q"[O2]*" g gDCW"^-1*" h"^-1),
  scales = list(alternating = FALSE),
  panel=function(x, y, ...) {
    panel.grid(h = -1, v = -1, col = grey(0.9))
    panel.errbars(x, y, ewidth = 0, ...)
    panel.key(..., which.panel = 2, corner = c(0.9, 0.9))
  }
)

print(plot_CO2_em, split = c(1,1,1,2), more = TRUE)
print(plot_O2_up, split = c(1,2,1,2))

3.7 High throughput imaging of chemostat cultivated cells

The same samples as were used to determine bulk PHB content using nile red fluorescence spectroscopy were also subjected to high throughput imaging. A Nikon microscope was used to acquire up to 5-10 images per replicate and condition. Phase contrast images and Texas Red fluorescence were acquired to identify cell outline/shape and Nile red fluorescnce from PHB staining, respectively. A custom Cell Profiler pipeline was used to automatically identfy cells, filter out non-cell objects by shape/size, and quantify fluorescence.

First we load a processed data frame combining data of fructose and ammonium limitation.

df_imaging <- read_csv("../data/input/20190728_chemostat_imaging.csv", col_types = cols())

Now we can plot the length of the cells and the volume calculated from the major axis length, and from the assumption that cells are rod/cylinder shaped objects. More precisely, cell volume using the volume of a cylinder as model was determined using the height h, (longest/major axis of our cell image), and its radius r (half of the minor axis/diameter of our cells). Then, volume is calculated as V = pi x r^2 x h.

Result: there does not seem to be a strong correlation of growth rate and cell size, neither in length nor volume. This is somewhat unexpected as this relationship was established in E. coli, Synechocystis, P. putida, and many other bacteria. However, the data might be incomplete and this question was not further pursued.

# see graph pars "box.rectangle" for the box, "box.umbrella" for whiskers
custom_boxplot <- custom.colorblind()
custom_boxplot$box.rectangle$lwd = 2
custom_boxplot$box.umbrella$lwd = 2

cell_length <- xyplot(MajorAxisLength_um ~ factor(Metadata_mu) | limitation, 
  df_imaging %>% filter(strain == "WT (PHB+)"),
  par.settings = custom_boxplot,
  scales = list(alternating = FALSE),
  xlab = expression("µ [h"^-1*"]"),
  ylab = expression("major axis length [µm]"),
  ylim = c(-0.5, 6.5),
  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.bwplot(x, y, horizontal = FALSE, pch = "|", 
      notch = FALSE, do.out = FALSE, ...)
    n_cells = tapply(y, x, length)
    panel.text(labels = paste("n > ", min(n_cells, na.rm = TRUE)-1), x = 2, y = 5.5)
  }
)

cell_volume <- xyplot(cell_volume_um3 ~ factor(Metadata_mu) | limitation, 
  df_imaging %>% filter(strain == "WT (PHB+)"),
  par.settings = custom_boxplot,
  scales = list(alternating = FALSE),
  xlab = expression("µ [h"^-1*"]"),
  ylab = expression("cell volume [µm"^3*"]"),
  ylim = c(-3, 33),
  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.bwplot(x, y, horizontal = FALSE, pch = "|", 
      notch = FALSE, do.out = FALSE, ...)
    n_cells = tapply(y, x, length)
    panel.text(labels = paste("n > ", min(n_cells, na.rm = TRUE)-1), x = 2, y = 28)
  }
)

print(cell_length, split = c(1,1,1,2), more = TRUE)
print(cell_volume, split = c(1,2,1,2))