03-model-analysis.qmd

---
filters: 
    - lightbox
lightbox: auto
---

# Models to the rescue; filamentation abstraction {#sec-model-analysis}

```{r}
#| label: libraries
#|
library(tidyverse)
library(patchwork)
library(here)
library(fs)
```

```{r}
#| label: set-default-plot-style
theme_set(
  theme_bw() +
  theme(
    legend.position = "top",
    strip.background = element_blank()
  )
)
```

## Introduction

By integrating information from the environment, cells can alter their
cell cycle. For instance, some cells arrest the cell division in the
presence of toxic agents but continue to grow. Previous studies have
shown that this filamentation phenomenon provides a mechanism that
enables cells to cope with stress, which increases the probability of
survival [@justiceMorphologicalPlasticityBacterial2008]. For example,
filamentation can be a process capable of subverting innate defenses
during urinary tract infection, facilitating the transition of
additional rounds of intracellular bacterial community formation
[@justiceFilamentationEscherichiaColi2006].

Although filament growth can help mitigate environmental stress (e.g.,
by activating the SOS response system
[@justiceMorphologicalPlasticityBacterial2008]), the evolutionary
benefits of producing elongated cells that do not divide are unclear.
Here, we proposed a mathematical model based on ordinary differential
equations that explicitly considers the concentration of intracellular
toxin as a function of the cell's length (see @eq-model-equation. The
model is built based on the growth ratio of measurements of the surface
area ($SA$) and the cell volume ($V$), whereby the uptake rate of the
toxin depends on the $SA$. However, $V's$ rate of change for $SA$ is
higher than $SA$ for $V$, which results in a transient reduction in the
intracellular toxin concentration (see
@fig-cell-dimensions-relationship)). Therefore, we hypothesized that
this geometric interpretation of filamentation represents a biophysical
defense line to increase the probability of a bacterial population's
survival in response to stressful environments.

```{r}
#| label: fig-cell-dimensions-relationship
#| fig-scap: Cell dimensions relationship.
#| fig-cap: >
#|   **Cell dimensions relationship.**
#|   We evaluated the resulting geometric properties on a grid of side lengths
#|   and radii with a pill-shaped cell. By maintaining a
#|   constant radius (typical case in bacteria such as *E. coli*) and
#|   increasing the side length, the surface area / volume relationship
#|   ($SA/V$) tends to decline since the $V$ will grow at a higher rate
#|   than the $SA$.
#|   
if (knitr::is_html_output()) {
  knitr::include_graphics("https://raw.githubusercontent.com/jvelezmagic/CellFilamentation/main/plots/cell_dimensions_relationship.svg")
} else {
  linguisticsdown::include_graphics2("https://raw.githubusercontent.com/jvelezmagic/CellFilamentation/main/plots/cell_dimensions_relationship.pdf")
}
```

## Filamentation model

Let us assume the shape of cells is a cylinder with hemispherical ends.
Based on this geometric structure, a nonlinear system of differential
equations governing filamentation can be written as follows:

$$
\begin{split}
\frac{dT_{int}}{dt} &= T_{sa} \cdot (T_{ext}(t) - T_{vol}) - \alpha \cdot T_{ant} \cdot T_{int} \\
\frac{dL}{dt} &=   \begin{cases}     \beta \cdot L,& \text{if } T_{int} \geq T_{sos},  t \geq \tau_{sos} + \tau_{delay} \text{ and } L < L_{max}  \\    0,            & \text{otherwise}  \end{cases}
\end{split}
$$ {#eq-model-equation}

It considers the internal toxin ($T_{int}$) and the cell length ($L$) as
variables. $T_{sa}$ and $T_{vol}$ represent the surface area and volume
of the toxin in the cell, respectively. $T_{ext}(t)$ is a function that
returns the amount of toxin in the cell medium. $T_{anti}$ and $\alpha$
symbolize the amount of antitoxin and its efficiency rate, respectively.
$\beta$ as the rate of filamentation. $L_{max}$ is the maximum size the
cell can reach when filamentation is on. $T_{sos}$ and $T_{kill}$ are
thresholds for filamentation and death, respectively. Finally,
$\tau_{delay}$ is the amount of time required to activate filamentation
after reaching the $T_{sos}$ threshold.

## Numerical results

### Filamentation provides transient resistance to stressful conditions

When growing rod-shaped bacterial cells under constant conditions, the
distribution of lengths and radii is narrow
[@schaechterGrowthCellNuclear1962]. However, some cells produce
filaments when growing under stress conditions
[@schaechterDependencyMediumTemperature1958]. Among the stress
conditions that can trigger the SOS response is exposure to beta-lactam
antibiotics [@millerSOSResponseInduction2004]. This phenomenon may
depend on the SOS response system
[@bosEmergenceAntibioticResistance2015], which can repair DNA damage,
giving the cell greater chances of recovering and surviving under stress
conditions. Besides, the SOS response has been reported to have precise
temporal coordination in individual bacteria
[@friedmanPreciseTemporalModulation2005].

In general, filamentation has been studied as an unavoidable consequence
of stress. However, we considered filamentation an active process that
produces the first line of defense against toxic agents. Therefore, a
differential equation model was proposed that assesses the change in the
amount of internal toxin as a function of cell length. At the core of
this model, we include the intrinsic relationship between the surface
area and the capsule volume since it is vital in determining cell size
[@harrisRelativeRatesSurface2016].

In @fig-filamentation-model-ramp-signal, cells grow in a ramp-shaped
external toxin signal without considering a toxin-antitoxin system. As
time progresses, the toxin in the external environment increases, so the
cell raises its internal toxin levels. At approximately time $22$, any
cell reaches $T_{sos}$. The control cell (unable to filament) and the
average cell (capable of filamenting) reach the death threshold,
$T_{kill}$, at times $31$ and $93$ (hatched and solid black lines),
respectively. Therefore, under this example, the cell has increased its
life span three times more than the control by growing as a filament
(green shaded area versus orange shaded area). In turn,
@fig-filamentation-model-ramp-signal also shows stochastic simulations
of the system in the intake of internal toxins. Since cell growth and
death processes are inherently stochastic, stochastic simulations would
be a better approximation. However, from now on, we will continue
studying the deterministic model.

```{r}
#| label: fig-filamentation-model-ramp-signal
#| fig-scap: Effect of filamentation on intracellular toxin concentration.
#| fig-cap: >
#|   **Effect of filamentation on intracellular toxin concentration.**
#|   In the presence of an extracellular toxic agent, the intracellular
#|   concentration of the toxin ($T_{int}$) increases until reaching a damage
#|   threshold that triggers filamentation ($T_{sos}$, blue point), increasing 
#|   cell length ($L$). When filamentation is on, the cell decreases $T_{int}$
#|   due to the intrinsic relationship between surface area and cell volume.
#|   When the cell reaches its maximum length, it dies if the
#|   stressful stimulus is not removed ($T_{kill}$, red dot). The hatched line
#|   represents a cell that can not grow as a filament.
#|   The orange shaded area represents the time between stress and the
#|   non-filament cell’s death, while the green shaded area represents the
#|   temporal gain. The blue background lines represent stochastic simulations
#|   of the same system.
#|
if(knitr::is_html_output()){
  knitr::include_graphics("https://raw.githubusercontent.com/jvelezmagic/CellFilamentation/main/plots/filamentation_model_ramp_signal.svg")
} else {
  linguisticsdown::include_graphics2("https://github.com/jvelezmagic/CellFilamentation/raw/main/plots/filamentation_model_ramp_signal.pdf") 
}
```

### Filamentation increases the minimum inhibitory concentration

```{r model-02-load-data}
#| results: hide
#| 
base_url <- "https://raw.githubusercontent.com/jvelezmagic/CellFilamentation/main/data/exp_raw/"
model_datasets <- c(
  "antitoxin_experiment.csv",
  "antitoxins_distributions.csv",
  "df_antitoxin_experiment.csv",
  "increase_resistance.csv",
  "toxin_exposure_experiment.csv"
) %>% 
  paste0(base_url, .) %>% 
  set_names(
    x = .,
    nm = . %>%
      path_file() %>% 
      path_ext_remove()
  ) %>% 
  map(read_csv, show_col_types = FALSE)

names(model_datasets)
```

In other to characterize the degree of resistance, dose-response
experiments determine the Minimum Inhibitory Concentration (MIC)
[@andrews2001; @andrewsDeterminationMinimumInhibitory2002]. Bacteria are
capable of modifying their MIC through various mechanisms, for example,
mutations [@lambertBacterialResistanceAntibiotics2005], impaired
membrane permeability [@satoOuterMembranePermeability1991], flux pumps
[@webberImportanceEffluxPumps2003], toxin-inactivating enzymes
[@wrightBacterialResistanceAntibiotics2005], and plasticity phenotypic
[@justiceMorphologicalPlasticityBacterial2008]. The latter is our
phenomenon of interest because it considers the change in shape and
size, allowing us to study it as a strategy to promote bacterial
survival.

We analyzed the MIC change caused by filamentation through stable
exposure experiments of different toxin amounts at other exposure times.
Computational simulations show that when comparing cells unable to
filament with those that can, there is an increase in the capacity to
tolerate more generous amounts of toxin, increasing their MIC (see
@fig-toxin-exposure-experiment). Therefore, it confers a gradual
increase in resistance beyond filamentation's role in improving the
cell's life span as the exposure time is longer.

```{r}
#| label: fig-toxin-exposure-experiment
#| fig-scap: Effect of filamentation on minimum inhibitory concentration (MIC).
#| fig-cap: >
#|   **Effect of filamentation on minimum inhibitory concentration (MIC).**
#|   By exposing a cell to different toxin concentrations with stable signals,
#|   the cell achieves a set MIC for conditions without or with filamentation
#|   (separation between stressed and dead state) for each exposure time,
#|   without representing a change for the normal state cells' points (blue zone). Thus,
#|   the green line represents a gradual MIC increase when comparing each MIC
#|   between conditions for each exposure time.
p_1 <- model_datasets$toxin_exposure_experiment %>%
  mutate(
    state = factor(state, levels = c("Normal", "Stressed", "Dead")),
    experiment = factor(experiment, levels = c("Control", "Normal"), labels = c("Without filamentation", "With filamentation"))
  ) %>% 
  ggplot(aes(x = exposure_time, y = amount_toxin, fill = state)) +
  geom_tile() +
  facet_grid(. ~ experiment) +
  scale_x_continuous(expand = c(0, 0)) +
  scale_y_continuous(expand = c(0, 0)) +
  scale_fill_hue(direction = -1, h.start = 90) +
  theme_minimal() +
  theme(
    legend.position = "top",
    panel.spacing.x = unit(1, "lines"),
    plot.tag.position = "topright"
  ) +
  labs(
    x = "",
    y = "Amount toxin",
    fill = "Cell status",
    tag = "A"
  ) +
  scale_fill_manual(
      values = c("#43b284", "#fab255", "#dd5129")
  ) +
  NULL

p_2 <- model_datasets$increase_resistance %>% 
  ggplot(aes(x = exposure_time, y = tolerance)) +
  geom_line() +
  geom_area(alpha = 1/ 3) +
  labs(
    x = "Exposure time",
    y = "Increased\nresistance",
    tag = "B"
  ) +
  scale_x_continuous(expand = c(0, 0)) +
  scale_y_continuous(expand = c(0, 0)) +
  theme_classic() +
  theme(
    plot.tag.position = "topright"
  )

p_toxin_exposure_experiment <- (p_1 / p_2) +
  plot_layout(heights = c(8,4))

p_toxin_exposure_experiment
```

### Heterogeneity in the toxin-antitoxin system represents a double-edged sword in survival probability

One of the SOS response system properties is that it presents
synchronous activation times within homogeneous populations
[@friedmanPreciseTemporalModulation2005]. It has constant gene
expression rates that help it cope with stress; however, it is possible
to introduce variability by considering having multicopy resistance
plasmids [@million-weaverMechanismsPlasmidSegregation2014]. Therefore,
the response times would have an asynchronous behavior at the global
level but synchronous at the local level.

To include this observation in the model, we include a negative term for
the internal toxin representing a toxin-antitoxin system. Therefore, the
model now has an efficiency rate of the antitoxin and a stable amount
per cell. We simulate the effect of the toxin-antitoxin system variation
within a $1000$-cell population; we initialize each one with similar
initial conditions, except for the amount of internal antitoxin, defined
as $T_{anti} \sim N(\mu, \sigma)$. Considering that $T_{anti}$ values
$< 0$ are equal to $0$. For each experiment, $\mu = 25$, while it was
evaluated in the range $[0-20]$. For the generation of pseudo-random
numbers and to ensure the results' reproducibility, the number $42$ was
considered seed.

As shown in @fig-antitoxin-experiment), when we compare heterogeneous
populations in their toxin-antitoxin system, we can achieve different
population dynamics, that is, changes in the final proportions of cell
states; normal, stressed, and dead. This difference is because the cell
sometimes has more or less antitoxin to handle the incoming stress.

```{r}
#| label: fig-antitoxin-experiment
#| fig-scap: Variability in the toxin-antitoxin system produces different proportions of cell states.
#| fig-cap: >
#|   **Variability in the toxin-antitoxin system produces different proportions of cell states.** 
#|   Histograms represent the distribution of antitoxin quantity, while the
#|   curves represent the population's fraction over time. The cell will start
#|   to filament after reaching a certain internal toxin threshold, $T_{sos}$.
#|   Therefore, the expected global effect on the population's response times
#|   based on the amount of antitoxin is asynchronous, while at the local
#|   level, it is synchronous. Consequently, different proportions are
#|   presented in the cellular states since some cells will activate the
#|   filamentation system before and others later.
#|
antitoxin_dist_p <- model_datasets$antitoxins_distributions %>% 
  ggplot(aes(x = antitoxin, fill = ..x..)) +
  geom_histogram(binwidth = 5, breaks = seq(0, 70, 5)) +
  facet_grid(. ~ sigma) +
  scale_fill_viridis_c(option = "inferno", labels = scales::percent) +
  labs(
    x = "Amount of antitoxin",
    y = "Number of cells"
  ) +
  theme(
    legend.position = "none",
    panel.spacing.x = unit(1, "lines"),
    strip.background = element_blank(),
    plot.tag.position = "topright",
    panel.grid = element_blank()
  )

antitoxin_lines_p <- model_datasets$df_antitoxin_experiment %>% 
  mutate(
    variable = factor(variable, levels = c("Normal", "Stressed", "Dead")),
  ) %>% 
  ggplot(aes(x = time, y = value, fill = variable)) +
  geom_area() +
  facet_grid(. ~ sigma) +
  scale_x_continuous(expand = c(0, 0)) +
  scale_y_continuous(expand = c(0, 0)) +
  scale_fill_hue(direction = -1, h.start = 90) +
  labs(
    x = "Exposure time",
    y = "Fraction population",
    fill = "Cell status"
  ) +
  scale_fill_manual(
      values = c("#43b284", "#fab255", "#dd5129")
  ) +
  theme(
    panel.spacing.x = unit(1, "lines"),
    strip.background = element_blank(),
    plot.tag.position = "topright"
  )
  

p_variability_toxin_antitoxin <- (antitoxin_dist_p / antitoxin_lines_p) +
  plot_layout(guides = "collect")

p_variability_toxin_antitoxin
```

Considering that the toxin-antitoxin system's variability can modify the
proportions of final cell states, a question arose about heterogeneity
levels' global behavior. To answer this question, we evaluate the
probability of survival for each population, defined by its distribution
of antitoxin levels. In this way, we characterized the population
survival probability function into three essential points according to
its effect: negative, invariant, and positive (see
@fig-survival-probability). Furthermore, these points are relative to
the homogeneous population's death time in question ($\tau_{kill}$):
when $t < \tau_{kill}$ will represent a detrimental effect on survival,
$t = \tau_{kill}$ will be independent of variability, and
$t > \tau_{kill}$ will be a beneficial point for survival. Therefore, we
concluded that the effect of heterogeneity in a toxin-antitoxin system
represents a double-edged sword.

```{r}
#| label: fig-survival-probability
#| fig-scap: Effect of variability on the toxin-antitoxin system.
#| fig-cap: >
#|   **Effect of variability on the toxin-antitoxin system.**
#|   The color of the heatmap is representative of the fraction of living cells
#|   at exposure time. The white vertical line represents the death time of
#|   the homogeneous population ($\tau_{kill}$). At $t < \tau_{kill}$, it is
#|   shown that the fraction of survivors decreases as the variability in the
#|   population increases. When $t = \tau_{kill}$, the variability does not
#|   affect the fraction of survivors, but it represents a percentage
#|   improvement for the homogeneous population. Finally, when
#|   $t > \tau_{kill}$, the heterogeneity of the population favors survival.
#|
survival_probability <- model_datasets$antitoxin_experiment %>% 
  mutate(
    value = value / max(value)
  ) %>% 
  filter(variable == "Dead") %>% 
  mutate(value = 1 - value) %>% 
  identity()
  
population_t_kill = survival_probability %>% 
  filter(sigma == 0, value == 0) %>% 
  pull(time) %>% 
  first()

p_survival_probability <- survival_probability %>% 
  ggplot(aes(x = time, y = sigma, fill = value)) +
  geom_raster() +
  geom_vline(xintercept = population_t_kill, color = "white", linetype = "dashed") +
  scale_x_continuous(expand = c(0, 0)) +
  scale_y_continuous(expand = c(0, 0)) +
  scale_fill_viridis_c(option = "inferno", labels = scales::percent) +
  labs(
    x = "Exposure time",
    y = "Population variability",
    fill = "Survival probability"
  ) +
  theme(
    legend.position = "top"
  ) +
  NULL

p_survival_probability
```

## Benefits and limitations of the model

Today, there have been advancements in the number of techniques that
have allowed it to extend quantitative analyses to individual cells'
dynamic observations [@camposConstantSizeExtension2014;
@meldrumFacultyOpinionsRecommendation2005;
@sliusarenkoHighthroughputSubpixelPrecision2011;
@taheri-araghiCellSizeControlHomeostasis2017;
@ursellRapidPreciseQuantification2017]. Therefore, studying their
cellular behavior daily from medium to medium can be somewhat
reproducible, facilitating the association of complex biological
functions in simple mathematical equations
[@neidhardtBacterialGrowthConstant1999].

Here, we proposed a mathematical model showing that filamentation could
be a population's resilience mechanism to stress conditions. Finding
that filamentation's net effect generates an additional window of time
for the cell to survive, decreasing the toxin's intracellular
concentration. However, we also found that filamentation's side effect
increases the cell's minimum inhibitory concentration. On the other
hand, when we introduce variability in the amount of antitoxin in a cell
population, we found that heterogeneity can be a double-edged sword,
sometimes detrimental and sometimes beneficial, depending on the time of
exposure to the toxic agent.

Notwithstanding the lack of parameters that are a little closer to
reality, confirming that the model can work under experimental
conditions would represent an achievement due to its explanatory
simplicity. Due to the above, despite being simple, the model could be
able to recapitulate the behavior seen in nature from variables that we
can easily calculate with single-cell measurements. However, in other
situations, it could be helpful to consider adding variables such as
cell wall production and peptidoglycans' accumulation, among others.
Starting in this way, the study of filamentation as a mechanism oriented
to the ecology of stress.