```{r include = FALSE, echo = FALSE}
knitr::opts_chunk$set(include = TRUE,
                      warning = FALSE,
                      message = FALSE,
                      cache = FALSE)
```

```{r echo = FALSE}
rm(list = ls())
library(tidyverse)
library(microxanox)
library(patchwork)
library(here)
source(here("r", "various_useful_functions.R"))
# source(here("experiments", "1_figures_main_ms", "ms_figure_functions.R"))
source(here("R/ms_figure_functions.R"))

#print(fig_div_vs_o2diff_1strain_7row)

colfunc_CB <- colorRampPalette(c("#024F17", "#B5FFC9"))
colfunc_SB <- colorRampPalette(c("#7D1402", "#FCBEB3"))
colfunc_PB <- colorRampPalette(c("#6E0172", "#F9AEFC"))
zero <- 0
unity <- 1
```

```{r}
microxanox_version_here <- "0.9.1"
microxanox_version_path <- paste0("microxanox_v", microxanox_version_here)
if (packageVersion("microxanox") < package_version(microxanox_version_here)) {
  stop(paste0("microxanox version needs to be at least ", microxanox_version_here))
}
print(paste0("Version of microxanox: ", microxanox_version_here))

```


# Introduction

This document is the Supplementary Report for the report *Functional diversity can facilitate the collapse of an undesirable ecosystem state*.


# Growth and inhibition functions

The rate of change in population densities is given by the per capita growth rate - per capita death rate, multiplied by population size ([Bush et al 2017](https://www.nature.com/articles/s41467-017-00912-x)). Growth rate is the growth rate in absence of any inhibition multiplied by the amount of inhibition. 

Growth rate in the absence of inhibition was modeled with the Monod equation with multiplicative kinetics when two substrates determine growth rate:

$g_j(X, Y) = g_{max,j} (\frac{X}{K_{j,X} + X})(\frac{Y}{K_{j,Y} + Y})$

$g_{max,j}$ is the maximum specific growth rate of strain $j$. $X$ is the concentration of substrate 1, and $Y$ is the concentration of substrate 2. $K_{j,X}$ and $K_{j,Y}$ are the half-saturation constants of species $j$ on substrates $X$ and $Y$.

Inhibition of growth by a substrate was modeled with the Haldane equation:

$h_j(X) = \frac{1}{1+ (X / H_{j,X})}$

$H_{j,X}$ can be thought of as a half-inhibition constant: the concentration of $X$ at which the growth rate of species $j$ is reduced by 50%. A higher half-inhibition constant therefore means higher tolerance of species $j$ to substrate $X$.

For further information on the functions used to model the system, please see the original publication by [Bush et al 2017](https://www.nature.com/articles/s41467-017-00912-x).

# Creating among strain diversity

Here we provide additional description of the method used to create strain variation within functional groups. The following should be read only after the relevant text in the main report.  We generated variation in maximum growth rate $g_{max}$ and in tolerance to reduced sulfur or oxygen (i.e. in the half-inhibition constants $H_{SR}$ or $H_{O}$) by varying the range of trait values around a reference trait value. To manipulate the amount of variation in a given trait and functional group, we introduced the meta-parameter $Div$ subscripted by the trait and the functional group. Take for example the cyanobacteria functional group, a reference maximum growth rate parameter $g_{max,CB}$ and a reference tolerance parameter $H_{SR,CB}$. The amount of among strain variation in the cyanobacteria is controlled by the two meta-parameters $Div_{g_{max,CB}}$ and $Div_{H_{SR,CB}}$. To create a trade-off, the two meta-parameters were always of different sign. The among strain variation was calculated such that the growth rate of strain $i = 1$ was a factor $1/({2Div_{g_{max,CB}}})$ of the reference growth rate, and the growth rate of the $i = n$ strain was a factor $2Div_{g_{max,CB}}$ of the reference growth rate. The growth rate of the $i = {2,...n-1}$ other strains was the reference growth rate multiplied by a factor that was equally distributed between the factor of the $i = 1$ and $i = n$ strain. This implementation of the variation and trade-off has the assumption of a linear trade-off of log-log transformed among-strain variation. Hence, for example, with the reference growth rate $g_{max,CB} = 0.05$, among strain variation of $Div_{g_{max,CB}} = 1$, and nine strains ($n = 9$), the growth rates of the nine strains are 0.025, 0.03, 0.035, 0.042, 0.05, 0.059, 0.071, 0.084, 0.1. Put another way, the growth rate parameter of the $i$-th of the $n$ strains in a functional group was calculated as $g_{max,i} = D_i * g_{max}$ where $D_i$ is the $i$-th element of the set $D = {2f(x)Div_{g_{max}} |x ∈ N, x ≤ n}$, where $f(x) = (x − (n + 1)/2)/((n − 1)/2)$.

This procedure ensures that the range of trait values is independent of the number of strains. When there were only two strains, the two strains took the minimum and maximum traits values, as defined by the $Div$ variable for that trait. When there were three strains, one had the minimum trait value, one the maximum, and one the mean. With nine strains, one had the minimum trait value, one the maximum, one the mean, and three were equally spaced (on log scale) between the minimum and the mean, and the remaining three between the maximum and the mean.

We show three examples of the created variation among strains within functional groups, one for nine strains (Figure \@ref(fig:nine-strain-div)), one for three (Figure \@ref(fig:three-strain-div)), and one for two (Figure \@ref(fig:two-strain-div)). In each of these three examples, three graphs show the diversity among strains in the maximum growth rate trait and the tolerance trait, plotted against each other and therefore also showing the trade-off between the two traits.

Notice that in all three cases, the range of trait variation represented by the strains is the same, even though the number of strains differs. This approach allowed us to independently manipulate trait diversity and number of strains. Also note that maximum growth rate of strains is evenly distributed on a log scale, even though it may appear at first site that they are distributed evenly on a linear scale (e.g. Figure \@ref(fig:three-strain-div)).

```{r  echo = FALSE}
CB_var_multiplier <- 2
SBPB_var_multiplier <- 6
CB_gmax_div <- 0.015789474 * CB_var_multiplier
CB_h_div <- -0.08 * CB_var_multiplier
SB_gmax_div <- 0.015789474 * SBPB_var_multiplier
SB_h_div <- -0.323  * SBPB_var_multiplier
PB_gmax_div <- 0.015789474  * SBPB_var_multiplier
PB_h_div <- -0.323  * SBPB_var_multiplier
num_div_treatment_levels <- 2
```

All three examples have the following amounts of trait variation: $Div_{g_{max,CB}}$ = `r round(CB_gmax_div, 3)`, $Div_{H_{SR,CB}}$ = `r round(CB_h_div, 3)`, $Div_{g_{max,SB}}$ = `r round(SB_gmax_div, 3)`, $Div_{H_{O,SB}}$ = `r round(SB_h_div, 3)`, $Div_{g_{max,PB}}$ = `r round(PB_gmax_div, 3)`, $Div_{H_{O,PB}}$ = `r round(PB_h_div, 3)`. These are also the maximum amounts of trait variation explored in the simulations.

(ref:nine-strain-div) Trait variation with nine strains in each of the functional groups. In all panels colour indicates the functional group, darker shades show more tolerant-lower growth rate strains, lighter shades show less tolerant-higher growth rate strains. Tolerance is given as the concentration of the inhibiting substrate where the growth rate of the strain is reduced by 50%. 

```{r nine-strain-div, echo = FALSE,  include=TRUE, fig.cap='(ref:nine-strain-div)', fig.align="center", fig.height=3, fig.width=8}

num_CB_strains <- 9
num_PB_strains <- 9
num_SB_strains <- 9
var_expt <- create_diversity_factorial(
  zero = zero, unity = unity,
  num_div_treatment_levels = num_div_treatment_levels,
  CB_gmax_div = CB_gmax_div, CB_h_div = CB_h_div,
  SB_gmax_div = SB_gmax_div, SB_h_div = SB_h_div,
  PB_gmax_div = PB_gmax_div, PB_h_div = PB_h_div,
  default_9strain = new_strain_parameter(
    n_CB = num_CB_strains,
    n_SB = num_SB_strains,
    n_PB = num_PB_strains,
    values_initial_state =  "bush_ssfig3"
  )
)
display_diversity(4,
  var_expt = var_expt,
  num_CB_strains = num_CB_strains,
  num_SB_strains = num_SB_strains,
  num_PB_strains = num_PB_strains,
  which_graphs = 1
)
```


(ref:three-strain-div) Trait variation with three strains in each of the functional groups. Note that the three displayed strains are the 1st, 5th, and 9th strains from the nine strain example shown in Figure \@ref(fig:nine-strain-div). In all panels colour indicates the functional group, darker shades show more tolerant-lower growth rate strains, lighter shades show less tolerant-higher growth rate strains. Tolerance is given as the concentration of the inhibiting substrate where the growth rate of the strain is reduced by 50%. 

```{r three-strain-div, echo = FALSE,  include=TRUE, fig.cap='(ref:three-strain-div)', fig.height=3, fig.width=8, fig.align="center"}

num_CB_strains <- 3
num_PB_strains <- 3
num_SB_strains <- 3
var_expt <- create_diversity_factorial(
  zero = zero, unity = unity,
  num_div_treatment_levels = num_div_treatment_levels,
  CB_gmax_div = CB_gmax_div, CB_h_div = CB_h_div,
  SB_gmax_div = SB_gmax_div, SB_h_div = SB_h_div,
  PB_gmax_div = PB_gmax_div, PB_h_div = PB_h_div,
  default_9strain = new_strain_parameter(
    n_CB = num_CB_strains,
    n_SB = num_SB_strains,
    n_PB = num_PB_strains,
    values_initial_state =  "bush_ssfig3"
  )
)
display_diversity(4,
  var_expt = var_expt,
  num_CB_strains = num_CB_strains,
  num_SB_strains = num_SB_strains,
  num_PB_strains = num_PB_strains,
  which_graphs = 1
)


#0.05*2^(CB_gmax_div) - 0.05
#0.05 - 0.05*2^(-CB_gmax_div)


```



(ref:two-strain-div) Trait variation with two strains in each of the functional groups. Note that the two displayed strains are the 1st and 9th strains from the nine strain example shown in Figure \@ref(fig:nine-strain-div). In all panels colour indicates the functional group, darker shades show more tolerant-lower growth rate strains, lighter shades show less tolerant-higher growth rate strains. Tolerance is given as the concentration of the inhibiting substrate where the growth rate of the strain is reduced by 50%. 

```{r two-strain-div, echo = FALSE,  include=TRUE, fig.cap='(ref:two-strain-div)', fig.height=3, fig.width=8, fig.align="center"}

num_CB_strains <- 2
num_PB_strains <- 2
num_SB_strains <- 2
var_expt <- create_diversity_factorial(
  zero = zero, unity = unity,
  num_div_treatment_levels = num_div_treatment_levels,
  CB_gmax_div = CB_gmax_div, CB_h_div = CB_h_div,
  SB_gmax_div = SB_gmax_div, SB_h_div = SB_h_div,
  PB_gmax_div = PB_gmax_div, PB_h_div = PB_h_div,
  default_9strain = new_strain_parameter(
    n_CB = num_CB_strains,
    n_SB = num_SB_strains,
    n_PB = num_PB_strains,
    values_initial_state =  "bush_ssfig3"
  )
)
display_diversity(4,
  var_expt = var_expt,
  num_CB_strains = num_CB_strains,
  num_SB_strains = num_SB_strains,
  num_PB_strains = num_PB_strains,
  which_graphs = 1
)
```

A consequence of this variation and of the trade-offs is shown in  Figure \@ref(fig:growth-rates), which shows the relationship between realised growth rate and the concentration in the environment of the respective inhibiting substance, e.g. reduced sulfur for the cyanobacteria. In these graphs only the least and most tolerant strains are shown. The lines cross due to the trade-off. At high concentrations of the inhibiting substance the more tolerant strain has higher realised growth rate, while at lower concentrations it is the less tolerant strain that has the higher realised growth rate.

(ref:growth-rates) Response of realised growth rate to variation in concentration of the inhibiting substrate for two strains, namely the most and the least tolerant strain. As always colour indicates the functional group, darker shades show more tolerant-lower growth rate strains, lighter shades show less tolerant-higher growth rate strains.

```{r growth-rates, echo = FALSE,  include=TRUE, fig.cap='(ref:growth-rates)', fig.height=3, fig.width=8, fig.align="center"}
display_diversity(4,
  var_expt = var_expt,
  num_CB_strains = num_CB_strains,
  num_SB_strains = num_SB_strains,
  num_PB_strains = num_PB_strains,
  which_graphs = 2
)
```


We have confirmed that when there is no within functional group strain variation, the dynamics of the model are not affected by the number of strains in each of the functional groups. This is a necessary condition for such a model ([Moisset de Espanés et al. 2021](https://onlinelibrary.wiley.com/doi/10.1111/ele.13775)), and gives confidence that our methods and their implementation are correct and robust.


# Method for finding stable states

## General method

We now discuss two approaches for finding how the stable state of the system varies along the gradient of oxygen diffusivity when there exists the possibility of multiple stable states. In the main article we describe and report results of the method that we term the *Temporal method*.

Here we also describe and give results of what we term the *Replication method*. In this, we made separate simulations for each of many values of oxygen diffusivity, ran the simulations for long enough to ensure transients had passed, and recorded the state. We ran multiple simulations per value of oxygen diffusivity, and started them with different initial conditions (either low or high total cyanobacterial abundance), so as to check for presence of alternate stable states. This method was used by [Bush et al 2017](https://www.nature.com/articles/s41467-017-00912-x). We term it the replication method since it is achieved by having many independent replicates of the ecosystem, each with a different oxygen diffusivity, and also a separate replicate for each of the different chosen initial conditions.

As is shown in the figures below, the temporal and replication methods result in very similar regions of bistability.

### Replication method, no diversity


```{r repl-meth1, echo=FALSE, eval = TRUE}
repl_meth <- readRDS(here("data",
                          microxanox_version_path,
                          "replication_method",
                          "event_definition_2",
                          "ss_data_9strains_sim_length1e+06_event_definition_2.RDS"))
temp_meth <- readRDS(here("data",
                          microxanox_version_path,
                          "temporal_method",
                          "event_definition_2",
                          "ss_data_9strains_waittime1e+06_event_definition_2.RDS"))
```



(ref:repl-meth2) Stable states found with the **replication method**. The left column of panels refers to the stable state when cyanobacteria density is at first high (favouring the oxic state). The right column refers to when cyanobacteria density is at first low (favouring the anoxic state). Points (which sometimes are so close as to appear as thick lines) show the stable state at the end of each period of constant oxygen diffusivity. Thin lines are for visualisation only, and join the points. Grey arrows show the oxygen diffusivity at which the system flips, and the direction of the flip. The length of the grey arrows is without meaning.

```{r repl-meth2, echo = FALSE,  include=TRUE, fig.cap='(ref:repl-meth2)', fig.height=10, fig.width=8, fig.align="center"}
fig_state_vs_o2diff_sidebyside_alternative_repmeth(repl_meth[1,])
```

### Temporal method, no diversity

(ref:temp-meth2) Stable states found with the **temporal method**. The left column of panels refers to the stable state when oxygen diffusivity is at first high (favouring the oxic state). The right column refers to when oxygen diffusivity is at first low (favouring the anoxic state). Points (which sometimes are so close as to appear as thick lines) show the stable state at the end of each period of constant oxygen diffusivity. Thin lines are for visualisation only, and join the points. Grey arrows show the oxygen diffusivity at which the system flips, and the direction of the flip. The length of the grey arrows is without meaning.

```{r temp-meth2, echo = FALSE,  include=TRUE, fig.cap='(ref:temp-meth2)', fig.height=10, fig.width=8, fig.align="center"}
fig_state_vs_o2diff_sidebyside_alternative(temp_meth[1,])
```


### Both methods overlaid, no diversity

(ref:both-meth1) Stable states of the two methods overlaid, the temporal method with a dashed line, and the replication method with a solid line.

```{r both-meth1, echo = FALSE,  include=TRUE, fig.cap='(ref:both-meth1)', fig.height=8, fig.width=8, fig.align="center"}

repl_method_1 <- repl_meth$ss_res[[1]] %>%
     mutate(time = NA,
            direction = ifelse(initial_N_CB == 1, "up", "down")) %>%
  select(names(temp_meth$ssfind_result[[1]]))

ss_result1 <- repl_method_1
ss_result2 <- temp_meth$ssfind_result[[1]]

plot_ss_result6_new(repl_method_1,
                    temp_meth$ssfind_result[[1]],
                    xlims = c(-8, 0))
```


## Effect of diversity, replication method



### Overall effects of diversity

The replication method shows qualitatively the same effects of functional diversity (Figure \@ref(fig:div-repl-meth1)) as those revealed by the temporal method. Namely

-   Introducing functional diversity in only the cyanobacteria increases the resilience of the oxic state and slightly decreases the resilience of the anoxic state. 
-   Introducing functional diversity in only the sulfate-reducing bacteria increases the resilience of the anoxic state and reduces the resilience of the oxic state.
-   Introducing functional diversity in only the phototrophic sulfur bacteria reduces the resilience of the anoxic state and has no effect on the resilience of the oxic state. 
-   Introducing functional diversity in more than one functional group reduces, reverses, or does not change the effect of diversity in individual functional groups. For example, introducing functional diversity in all three functional groups reduces the resilience of the oxic state when trait variation is low but increases the resilience of the oxic state when trait variation is high. 

(ref:div-repl-meth1) **Replication method results.** **Left column of panels**: the effect of trait diversity (y-axis) on the positions of the tipping points (anoxic to oxic in blue, oxic to anoxic in red), which determine the extent of the region of bistability (green shaded region) in the oxygen diffusivity gradient. The left-most region is anoxic only, the middle is oxic or anoxic depending on historical conditions, and the right-most is oxic only. The thick grey vertical lines at -8 and 0 on the x-axis show the extent of the oxygen diffusivity gradient investigated in the simulations. In panels e, and i, the orange region shows when patterns of coexistence were atypical of the classical pattern shown in Figure 1a of the main text (e.g. when cyanobacteria coexisted with sulfur bacteria at low oxygen diffusivities). **Right column of panels**: The effect size of trait variation on resilience, with positive values on the y-axis indicating that diversity has increased resilience.

```{r div-repl-meth1, echo = FALSE,  include=TRUE, fig.cap='(ref:div-repl-meth1)', fig.height=8, fig.width=8, fig.align="center"}
all_stab <- readRDS(here("data", 
                         microxanox_version_path,
                         "replication_method",
                         "event_definition_2",
                         "processed_data",
                         "all_stab_data.RDS"))

sim_length <- 1e6 ## currently a problem in how sim length is parsed by the data processing step
num_strains <- 9

plot_here <- all_stab %>%
  mutate(stab_states_type = "standard",
         stab_states_type = case_when(
          (var_treat == "CB" & stand_var > 1.00) ~ "non-standard",
          (var_treat == "SB" & stand_var > 1.00) ~ "non-standard",
          (var_treat == "PB" & stand_var > 0.69) ~ "non-standard",
          (var_treat == "CB-SB" & stand_var > 1.00) ~ "non-standard",
          (var_treat == "CB-PB" & stand_var > 0.65) ~ "non-standard",
          (var_treat == "SB-PB" & stand_var > 1.00) ~ "non-standard",
          (var_treat == "CB-SB-PB" & stand_var > 1.00) ~ "non-standard"
         ),
         stab_states_type = ifelse(is.na(stab_states_type), "standard", "non-standard")
  )
       
  
  
p1 <- fig_div_vs_o2diff_1strain_7row_stable_state_type(plot_here,
                                     which_strain = num_strains,
                                     figure_title = paste0(sim_length,
                                                           " simulation length,\n",
                                                           num_strains,
                                                           " strains"))
p2 <- fig_resilience_vs_div(plot_here,
                            which_strain = num_strains,
                            figure_title = " ")
p1 + p2

```

### An example of stable states

Figure \@ref(fig:repl-meth-ex1) shows an example of stable states along the oxygen diffusivity gradient produced by the replication method for finding stable states. The patterns are similar to those produced by the temporal method for finding the stable states (see Figure 4 in the main report), though they differ in one respect: close to the transition back to the state a group dominates, some of the more tolerant strains occur in the stable states found by the temporal method (Figure 4 in the main report), while they do not occur in stable states found by the replication method (Figure \@ref(fig:repl-meth-ex1) here). We understand that the appearance of these strains is transient, and the difference in their appearance results from the differences in starting conditions of the two methods. Specifically, in the temporal method, the initial condition at a new level of oxygen diffusivity is the stable state at the previous level of oxygen diffusivity, while in the replication method the starting conditions are identical at each level of oxygen diffusivity. These differences in starting conditions apparently result in longer duration of transients in the temporal method than the replication method.

An interesting feature of these results is that although the position of the oxic to anoxic tipping point is affected by diversity in the sulfate-reducing bacteria (Figure \@ref(fig:div-repl-meth1) and Figure 3 in the main report) there is no evidence of more tolerant strains of sulfate-reducing bacteria being present in the final system state (Figure \@ref(fig:repl-meth-ex1)). This phenomenon is explained in Section \@ref(SB2).

(ref:repl-meth-ex1) Stable states found with the **replication method** with intermediate levels of diversity in all functional groups. Figure features are otherwise the same as those in Figure \@ref(fig:repl-meth2)).

```{r repl-meth-ex1, echo = FALSE,  include=TRUE, fig.cap='(ref:repl-meth-ex1)', fig.height=10, fig.width=8, fig.align="center"}
repl_meth <- readRDS(here("data",
                          microxanox_version_path,
                          "replication_method",
                          "event_definition_2",
                          "ss_data_9strains_sim_length1e+06_event_definition_2.RDS"))
fig_state_vs_o2diff_sidebyside_alternative(repl_meth[94,])
```


## Effect of abundance change event


As mentioned in the main report, to avoid very low strain densities, we added 1 unit of density to every biological state variable every 1,000 hours. We here show that a slightly different method produced the same results. Specifically, we tested the effect of making an abundance floor of 1 unit that was implemented every 1,000 hours. I.e. every 1,000 hours, any abundances that were lower than 1 unit were set to 1 unit. As shown in Figure \@ref(fig:div-temp-eventdef1), the results with this floor abundance are similar as when 1 was added every 1,000 hours. 

(ref:div-temp-eventdef1) **Results with an abundance floor.** **Left column of panels**: the effect of trait diversity (y-axis) on the positions of the tipping points (anoxic to oxic in blue, oxic to anoxic in red), which determine the extent of the region of bistability (green shaded region) in the oxygen diffusivity gradient. The left-most region is anoxic only, the middle is oxic or anoxic depending on historical conditions, and the right-most is oxic only. The thick grey vertical lines at -8 and 0 on the x-axis show the extent of the oxygen diffusivity gradient investigated in the simulations. In panels e, and i, the orange region shows when patterns of coexistence were atypical of the classical pattern shown in Figure 1a of the main text (e.g. when cyanobacteria coexisted with sulfur bacteria at low oxygen diffusivities). **Right column of panels**: The effect size of trait variation on resilience, with positive values on the y-axis indicating that diversity has increased resilience, i.e. delayed the shift to the other state.


```{r div-temp-eventdef1, eval = TRUE, echo = FALSE,  include=TRUE, fig.cap='(ref:div-temp-eventdef1)', fig.height=8, fig.width=8, fig.align="center"}
all_stab <- readRDS(here("data",
                          microxanox_version_path,
                          "temporal_method",
                          "event_definition_1",
                          "processed_data",
                         "all_stab_data.RDS"))

wait_time <- 1e6
num_strains <- 9

plot_here <- all_stab %>%
  mutate(stab_states_type = "standard",
         stab_states_type = case_when(
          (var_treat == "CB" & stand_var > 1.00) ~ "non-standard",
          (var_treat == "SB" & stand_var > 1.00) ~ "non-standard",
          (var_treat == "PB" & stand_var > 0.65) ~ "non-standard",
          (var_treat == "CB-SB" & stand_var > 1.00) ~ "non-standard",
          (var_treat == "CB-PB" & stand_var > 0.65) ~ "non-standard",
          (var_treat == "SB-PB" & stand_var > 1.00) ~ "non-standard",
          (var_treat == "CB-SB-PB" & stand_var > 1.00) ~ "non-standard"
         ),
         stab_states_type = ifelse(is.na(stab_states_type), "standard", "non-standard")
  )
       
  
  
p1 <- fig_div_vs_o2diff_1strain_7row_stable_state_type(plot_here,
                                     which_strain = num_strains,
                                     figure_title = paste0(sim_length,
                                                           " simulation length,\n",
                                                           num_strains,
                                                           " strains"))
p2 <- fig_resilience_vs_div(plot_here,
                            which_strain = num_strains,
                            figure_title = " ")
p1 + p2






```



## Further thoughts

Here we give some further thoughts regarding the comparison of the temporal and replication method for stable state finding.

A particularly notable difference between the two approaches is in the initial condition. In the temporal approach, the initial condition at a new level of oxygen diffusivity is the stable state at the previous level of oxygen diffusivity, while in the approach of [Bush et al 2017](https://www.nature.com/articles/s41467-017-00912-x) the initial conditions are fixed and arbitrary (though this approach is still useful for illustrating the bistability of the system).

The temporal approach is particularly useful when investigating the effects of strain richness. In preliminary simulations using the [Bush et al 2017](https://www.nature.com/articles/s41467-017-00912-x) approach, we found the counterintuitive result that higher strain richness led to lower resilience. We then recalled that we chose to divide the initial total abundance of the cyanobacterial (CB) functional group equally among the strains (termed a substitutive design in biodiversity manipulation experiments). This meant that the most tolerant CB strain in a three-strain simulation had higher initial abundance than the same most tolerant CB strain in a nine-strain simulation. The greater abundance of the most CB tolerant strain made it possible for a three-strain system to develop to be oxic, whereas the nine-strain system developed to be anoxic, despite identical range of trait variation. We could have chosen to do otherwise, for example start with the least tolerant strains of cyanobacteria as more abundant when the system started with high cyanobacterial abundance, and the most tolerant strain as more abundant when the system starts with low cyanobacterial abundance. Using the temporal approach removed the need to specify initial abundances (except the very first one), and therefore results were not so determined by the assumptions used to determine the initial conditions.


# Finding the separatrix

In the original publication ([Bush et al 2017](<https://www.nature.com/articles/s41467-017-00912-x>)) Figure 3a shows a separatrix in the cyanobacteria abundance - oxygen diffusivity dimensions of the system. If the system has initial cyanobacterial abundance greater than the value at the separatrix then the cyanobacterial abundance increases, and the system goes to the oxic stable state. Conversely, if the system has initial cyanobacterial abundance lower than the value at the separatrix then the cyanobacterial abundance decreases, and the system goes to the anoxic stable state.

It is critical to note, however, that the position of the separatrix depends on both the parameters of the system, and the initial conditions of variables other than cyanobacterial abundance. If the system were started with higher initial concentration, then the position of the separatrix would change (as would the values of oxygen diffusivity at which the system tipped from anoxic to oxic, and from oxic to anoxic).

This issue of dependence of the position of the separatrix on the value of initial condition of multiple state variables has particular significance when there are multiple strains per functional group. Then one cannot just vary initial total cyanobacterial abundance (as was done to create the separatrix on Figure 3a in the original publication) but one must also define how the initial abundance of each of the strains varies. One might choose to make each strain have equal initial abundance, and to vary the total (as we did in our implementation of the replication method in the main article).

We could, however, and arguably should, give the most tolerant strain of cyanobacteria higher initial relative abundance when total initial cyanobacterial abundance was low, and the least tolerant strain higher initial relative abundance when total cyanobacterial abundance was high. And we would need to specify the total and relative abundance of all cyanobacteria (and sulfur bactria) strains. Thus, with a system as complex as the 9 strain system, which has 31 state variables the separatrix is very high dimensional, and evaluating (and plotting) it in one dimension requires many assumptions. We chose to not attempt to find a meaningful set of assumptions in this study, though not because we think it would be unprofitable to do so, but rather because we did not want to risk distracting attention from findings we consider to be significant enough in their own right.

# Effect of simulation duration

## Diversity effects

Recall that the temporal method for finding stable states involves running the simulation for a certain amount of time, termed the "wait time", at a constant level of oxygen diffusivity. At the end of that wait time the state of the system is recorded, and then the oxygen diffusivity is slightly increased or decreased, and the system then simulated again for the duration of the wait time.

The following figures \@ref(fig:div-temp-meth1)-\@ref(fig:div-temp-meth4) show the effects of diversity with long to shorter wait times.

(ref:div-temp-meth1) Effects of diversity with wait time of 1e6 hours (same figure as in main report). (All else as in (Figure \@ref(fig:div-repl-meth1)).

```{r div-temp-meth1, echo = FALSE,  include=TRUE, fig.cap='(ref:div-temp-meth1)', fig.height=7, fig.width=6, fig.align="center"}
all_stab <- readRDS(here("data",
                         microxanox_version_path,
                         "temporal_method",
                         "event_definition_2",
                         "processed_data",
                         "all_stab_data.RDS"))

wait_time <- 1e6
num_strains <- 9


plot_here <- all_stab %>%
  filter(waittime == wait_time) %>%
  mutate(stab_states_type = "standard",
         stab_states_type = case_when(
          (var_treat == "CB" & stand_var > 1.00) ~ "non-standard",
          (var_treat == "SB" & stand_var > 1.00) ~ "non-standard",
          (var_treat == "PB" & stand_var > 0.65) ~ "non-standard",
          (var_treat == "CB-SB" & stand_var > 1.00) ~ "non-standard",
          (var_treat == "CB-PB" & stand_var > 0.65) ~ "non-standard",
          (var_treat == "SB-PB" & stand_var > 1.00) ~ "non-standard",
          (var_treat == "CB-SB-PB" & stand_var > 1.00) ~ "non-standard"
         ),
         stab_states_type = ifelse(is.na(stab_states_type), "standard", "non-standard")
  )
       
  
  
p1 <- fig_div_vs_o2diff_1strain_7row_stable_state_type(plot_here,
                                     which_strain = num_strains,
                                     figure_title = paste0(sim_length,
                                                           " simulation length,\n",
                                                           num_strains,
                                                           " strains"))
p2 <- fig_resilience_vs_div(plot_here,
                            which_strain = num_strains,
                            figure_title = " ")
p1 + p2

```

(ref:div-temp-meth2) Effects of diversity with wait time of 1e5 hours. (All else as in (Figure \@ref(fig:div-repl-meth1)).

```{r div-temp-meth2, echo = FALSE,  include=TRUE, fig.cap='(ref:div-temp-meth2)', fig.height=7, fig.width=6, fig.align="center"}

wait_time <- 1e5
num_strains <- 9


plot_here <- all_stab %>%
  filter(waittime == wait_time) %>%
  mutate(stab_states_type = "standard",
         stab_states_type = case_when(
          (var_treat == "CB" & stand_var > 1.00) ~ "non-standard",
          (var_treat == "SB" & stand_var > 1.00) ~ "non-standard",
          (var_treat == "PB" & stand_var > 0.65) ~ "non-standard",
          (var_treat == "CB-SB" & stand_var > 1.00) ~ "non-standard",
          (var_treat == "CB-PB" & stand_var > 0.65) ~ "non-standard",
          (var_treat == "SB-PB" & stand_var > 1.00) ~ "non-standard",
          (var_treat == "CB-SB-PB" & stand_var > 1.00) ~ "non-standard"
         ),
         stab_states_type = ifelse(is.na(stab_states_type), "standard", "non-standard")
  )
       
  
  
p1 <- fig_div_vs_o2diff_1strain_7row_stable_state_type(plot_here,
                                     which_strain = num_strains,
                                     figure_title = paste0(sim_length,
                                                           " simulation length,\n",
                                                           num_strains,
                                                           " strains"))
p2 <- fig_resilience_vs_div(plot_here,
                            which_strain = num_strains,
                            figure_title = " ")
p1 + p2

```

(ref:div-temp-meth3) Effects of diversity with wait time of 1e4 hours. (All else as in (Figure \@ref(fig:div-repl-meth1)).

```{r div-temp-meth3, echo = FALSE,  include=TRUE, fig.cap='(ref:div-temp-meth3)', fig.height=7, fig.width=6, fig.align="center"}

wait_time <- 1e4
num_strains <- 9


plot_here <- all_stab %>%
  filter(waittime == wait_time) %>%
  mutate(stab_states_type = "standard",
         stab_states_type = case_when(
          (var_treat == "CB" & stand_var > 1.00) ~ "non-standard",
          (var_treat == "SB" & stand_var > 1.00) ~ "non-standard",
          (var_treat == "PB" & stand_var > 0.65) ~ "non-standard",
          (var_treat == "CB-SB" & stand_var > 1.00) ~ "non-standard",
          (var_treat == "CB-PB" & stand_var > 0.65) ~ "non-standard",
          (var_treat == "SB-PB" & stand_var > 1.00) ~ "non-standard",
          (var_treat == "CB-SB-PB" & stand_var > 1.00) ~ "non-standard"
         ),
         stab_states_type = ifelse(is.na(stab_states_type), "standard", "non-standard")
  )
       
  
  
p1 <- fig_div_vs_o2diff_1strain_7row_stable_state_type(plot_here,
                                     which_strain = num_strains,
                                     figure_title = paste0(sim_length,
                                                           " simulation length,\n",
                                                           num_strains,
                                                           " strains"))
p2 <- fig_resilience_vs_div(plot_here,
                            which_strain = num_strains,
                            figure_title = " ")
p1 + p2

```

(ref:div-temp-meth4) Effects of diversity with wait time of 1e3 hours. (All else as in Figure \@ref(fig:div-repl-meth1)).

```{r div-temp-meth4, echo = FALSE,  include=TRUE, fig.cap='(ref:div-temp-meth4)', fig.height=7, fig.width=6, fig.align="center"}
wait_time <- 1e3
num_strains <- 9


plot_here <- all_stab %>%
  filter(waittime == wait_time) %>%
  mutate(stab_states_type = "standard",
         stab_states_type = case_when(
          (var_treat == "CB" & stand_var > 1.00) ~ "non-standard",
          (var_treat == "SB" & stand_var > 1.00) ~ "non-standard",
          (var_treat == "PB" & stand_var > 0.65) ~ "non-standard",
          (var_treat == "CB-SB" & stand_var > 1.00) ~ "non-standard",
          (var_treat == "CB-PB" & stand_var > 0.65) ~ "non-standard",
          (var_treat == "SB-PB" & stand_var > 1.00) ~ "non-standard",
          (var_treat == "CB-SB-PB" & stand_var > 1.00) ~ "non-standard"
         ),
         stab_states_type = ifelse(is.na(stab_states_type), "standard", "non-standard")
  )
       
  
  
p1 <- fig_div_vs_o2diff_1strain_7row_stable_state_type(plot_here,
                                     which_strain = num_strains,
                                     figure_title = paste0(sim_length,
                                                           " simulation length,\n",
                                                           num_strains,
                                                           " strains"))
p2 <- fig_resilience_vs_div(plot_here,
                            which_strain = num_strains,
                            figure_title = " ")
p1 + p2

```


```{r div-temp-meth5, echo = FALSE, eval = FALSE, include=TRUE, fig.cap='(ref:div-temp-meth5)', fig.height=7, fig.width=6, fig.align="center"}

## Owen removed these results from the report, because the tipping point finding algorithm was being misled
## by transient dynamics, and therefore was not finding the correct position of the anoxic to oxic tipping
## point

## (ref:div-temp-meth5) Effects of diversity with wait time of 1e2 hours. (All else as in Figure \@ref(fig:div-repl-meth1)).


wait_time <- 1e2
num_strains <- 9


plot_here <- all_stab %>%
  filter(waittime == wait_time) %>%
  mutate(stab_states_type = "standard",
         stab_states_type = case_when(
          (var_treat == "CB" & stand_var > 1.00) ~ "non-standard",
          (var_treat == "SB" & stand_var > 1.00) ~ "non-standard",
          (var_treat == "PB" & stand_var > 0.65) ~ "non-standard",
          (var_treat == "CB-SB" & stand_var > 1.00) ~ "non-standard",
          (var_treat == "CB-PB" & stand_var > 0.65) ~ "non-standard",
          (var_treat == "SB-PB" & stand_var > 1.00) ~ "non-standard",
          (var_treat == "CB-SB-PB" & stand_var > 1.00) ~ "non-standard"
         ),
         stab_states_type = ifelse(is.na(stab_states_type), "standard", "non-standard")
  )
       
  
  
p1 <- fig_div_vs_o2diff_1strain_7row_stable_state_type(plot_here,
                                     which_strain = num_strains,
                                     figure_title = paste0(sim_length,
                                                           " simulation length,\n",
                                                           num_strains,
                                                           " strains"))
p2 <- fig_resilience_vs_div(plot_here,
                            which_strain = num_strains,
                            figure_title = " ")
p1 + p2

```



## 1e6 vs 2e6 hours wait time

The results below show the stable states as a function of oxygen diffusivity for a system with an intermediate level of diversity in all functional groups, with either 1'000'000 hours at each oxygen diffusivity level (Figure \@ref(fig:ss-temp-meth3)), or 2'000'000 hours (Figure \@ref(fig:ss-temp-meth4)). This shows that 1'000'000 hours duration of the simulations are sufficient to reach stable state, except in regions of strain replacement, including at some state transitions. E.g. in the region of strain replacement in the sulfate-reducing bacteria there is difference in strain abundances at 1'000'000 and 2'000'000 hours. However, the position of the tipping points (which were calculated from the total SB abundance) and the substrate concentrations did not differ between 1'000'000 and 2'000'000 hours.

(ref:ss-temp-meth3) Stable states found with the temporal method with duration at each oxygen diffusivity value of **1e6 hours**. The left column of panels refers to the stable state when oxygen diffusivity is at first high (favouring the oxic state). The right column refers to when oxygen diffusivity is at first low (favouring the anoxic state). Points (which sometimes are so close as to appear as thick lines) show the stable state at the end of each period of constant oxygen diffusivity. Thin lines are for visualisation only, and join the points. Grey arrows show the oxygen diffusivity at which the system flips, and the direction of the flip. The length of the grey arrows is without meaning.


```{r ss-temp-meth3, eval = TRUE, echo = FALSE,  include=TRUE, fig.cap='(ref:ss-temp-meth3)', fig.height=10, fig.width=8, fig.align="center"}

temp_meth <- readRDS(here("data",
                          microxanox_version_path,
                          "temporal_method",
                          "event_definition_2",
                          "ss_data_9strains_waittime1e+06_event_definition_2.RDS"))

## data made with the "run_single_sim.R" file in the ss_find, temporal method folder.

single_1e6 <- readRDS(here("data",
                           microxanox_version_path, 
                           "temporal_method",
                           "single_sims",
                           "single_1e6_result.RDS"))
temp_meth[1,]$ssfind_result[[1]] <- single_1e6
fig_state_vs_o2diff_sidebyside_alternative(temp_meth[1,])
```


(ref:ss-temp-meth4) Stable states found with the temporal method with duration at each oxygen diffusivity value of **2e6 hours**. All other aspects are as in Figure \@ref(fig:ss-temp-meth3).


```{r ss-temp-meth4, eval = TRUE, echo = FALSE,  include=TRUE, fig.cap='(ref:ss-temp-meth4)', fig.height=10, fig.width=8, fig.align="center"}
single_2e6 <- readRDS(here("data",
                           microxanox_version_path, 
                           "temporal_method",
                           "single_sims",
                           "single_2e6_result.RDS"))
temp_meth[1,]$ssfind_result[[1]] <- single_2e6
fig_state_vs_o2diff_sidebyside_alternative(temp_meth[1,])
```

# Effect of number of strains

In the main text, we present results of simulations with nine strains per functional group. Here we also explore the effect of number of strains (2, 3, 6, 9) on the position of the two tipping points (Figure \@ref(fig:div-stab-multstrain1)). At low to moderate amounts of trait variation, the number of strains did not affect the position of the tipping points. However, at high amounts of trait variation, resilience was lower when the functional groups had only two or three strains rather than six or nine. This effect was particularly pronounced when there was high variation in only SB or in SB-CB: here, simulations with two or three strains led to dramatically lower resilience of the anoxic state than simulations with six or nine strains. We believe that this effect is due to the greater distance in trait space between strains in simulations with two or three strains than in simulations with more strains (also see Section 9). 

(ref:div-stab-multstrain1)  Effects of number of strains and of diversity among strains. In panels PB and CB-PB the region with dashed lines shows when patterns of coexistence were atypical of the classical pattern shown in Figure 1a of the main text (e.g. when cyanobacteria coexisted with sulfur bacteria at low oxygen diffusivities).

```{r div-stab-multstrain1, echo = FALSE,  include=TRUE, fig.cap='(ref:div-stab-multstrain1)', fig.height=8, fig.width=8, fig.align="center"}
all_stab <- readRDS(here("data",
                         microxanox_version_path,
                         "temporal_method", 
                         "event_definition_2",
                         "processed_data", "all_stab_data.RDS"))
which_strains <- c(2, 3, 6, 9)
which_wait_time<- 1e6
figure_title <- " "


plot_here111 <- all_stab %>%
  mutate(stab_states_type = "standard",
         stab_states_type = case_when(
          (var_treat == "CB" & stand_var > 1.00) ~ "non-standard",
          (var_treat == "SB" & stand_var > 1.00) ~ "non-standard",
          (var_treat == "PB" & stand_var > 0.65) ~ "non-standard",
          (var_treat == "CB-SB" & stand_var > 1.00) ~ "non-standard",
          (var_treat == "CB-PB" & stand_var > 0.65) ~ "non-standard",
          (var_treat == "SB-PB" & stand_var > 1.00) ~ "non-standard",
          (var_treat == "CB-SB-PB" & stand_var > 1.00) ~ "non-standard"
         ),
         stab_states_type = ifelse(is.na(stab_states_type), "standard", "non-standard")
  )
       

# plot_here111 <- all_stab %>%
#   filter(!(num_strains == 9 & var_treat == "CB" & stand_var > 1.00),
#          !(num_strains == 9 & var_treat == "SB" & stand_var > 1.00),
#          !(num_strains == 9 & var_treat == "PB" & stand_var > 0.65),
#          !(num_strains == 9 & var_treat == "CB-SB" & stand_var > 1.00),
#          !(num_strains == 9 & var_treat == "CB-PB" & stand_var > 0.65),
#          !(num_strains == 9 & var_treat == "SB-PB" & stand_var > 1.0),
#          !(num_strains == 9 & var_treat == "CB-SB-PB" & stand_var > 1.00),
#          
#          !(num_strains == 6 & var_treat == "CB" & stand_var > 1.00),
#          !(num_strains == 6 & var_treat == "SB" & stand_var > 1.00),
#          !(num_strains == 6 & var_treat == "PB" & stand_var > 0.65),
#          !(num_strains == 6 & var_treat == "CB-SB" & stand_var > 1.00),
#          !(num_strains == 6 & var_treat == "CB-PB" & stand_var > 0.65),
#          !(num_strains == 6 & var_treat == "SB-PB" & stand_var > 1.00),
#          !(num_strains == 6 & var_treat == "CB-SB-PB" & stand_var > 1.00),
#          
#          !(num_strains == 3 & var_treat == "CB" & stand_var > 1.00),
#          !(num_strains == 3 & var_treat == "SB" & stand_var > 1.00),
#          !(num_strains == 3 & var_treat == "PB" & stand_var > 0.65),
#          !(num_strains == 3 & var_treat == "CB-SB" & stand_var > 1.00),
#          !(num_strains == 3 & var_treat == "CB-PB" & stand_var > 0.65),
#          !(num_strains == 3 & var_treat == "SB-PB" & stand_var > 1.0),
#          !(num_strains == 3 & var_treat == "CB-SB-PB" & stand_var > 1.00),
#          
#          !(num_strains == 2 & var_treat == "CB" & stand_var > 1.00),
#          !(num_strains == 2 & var_treat == "SB" & stand_var > 1.00),
#          !(num_strains == 2 & var_treat == "PB" & stand_var > 0.65),
#          !(num_strains == 2 & var_treat == "CB-SB" & stand_var > 1.00),
#          !(num_strains == 2 & var_treat == "CB-PB" & stand_var > 0.65),
#          !(num_strains == 2 & var_treat == "SB-PB" & stand_var > 1.00),
#          !(num_strains == 2 & var_treat == "CB-SB-PB" & stand_var > 1.00))

fig_div_vs_o2diff_multistrain_stable_state_type(plot_here111, which_strains,
                              which_wait_time, figure_title)
```



# CB diversity effects

Increased CB diversity causes a slight decrease in the resilience of the anoxic state. The following two figures show the state responses to oxygen diffusivity with low (Figure \@ref(fig:ss-CB-only1)) and high (Figure \@ref(fig:ss-CB-only2)) levels of CB diversity. The difference in the position of the anoxic to oxic tipping point is small. It is, in fact the smallest it could be, as it is the smallest step size in oxygen diffusivity used in the simulations. 

(ref:ss-CB-only1) Stable state of the system with **low** level of diversity in only CB, and no diversity in the other two groups. All other aspects are as in Figure \@ref(fig:ss-temp-meth3).

```{r ss-CB-only1, echo = FALSE,  include=TRUE, fig.cap='(ref:ss-CB-only1)', fig.height=10, fig.width=8, fig.align="center"}
wait_time <- "1e+06"
ss_9s <- readRDS(here::here("data",
                            microxanox_version_path, 
                            "temporal_method",
                            "event_definition_2",
                            paste0("ss_data_9strains_waittime", wait_time, "_event_definition_2.RDS")))
#stab_9s <- readRDS(here::here(\"data", "0_ss_finding", "temporal_method", microxanox_version_path, paste0("stab_data_9strains_waittime", wait_time, "_event_definition_2.RDS"))


#max <- max(ss_3s$PB_var_gmax_s)
#f <- 0.625
#fv <- f*max

only_CB_variable <- ss_9s %>%
  filter(
    abs(SB_var_gmax_s) < 0.0001,
    abs(PB_var_gmax_s) < 0.0001
  )
only_CB_variable <- only_CB_variable[order(only_CB_variable$CB_var_gmax_s), ]

#only_CB_variable$CB_var_gmax_s

fig_state_vs_o2diff_sidebyside_alternative(only_CB_variable[2,])

#fig_state_vs_o2diff(only_SB_variable[13,])
```

(ref:ss-CB-only2) Stable state of the system with **high** level of diversity in only CB, and no diversity in the other two groups. All other aspects are as in Figure \@ref(fig:ss-temp-meth3).

```{r ss-CB-only2, echo = FALSE,  include=TRUE, fig.cap='(ref:ss-CB-only2)', fig.height=10, fig.width=8, fig.align="center"}

fig_state_vs_o2diff_sidebyside_alternative_leftfliponly(only_CB_variable[20,])

#ss_result <- only_CB_variable[20,]

#fig_state_vs_o2diff(only_SB_variable[13,])
```



# SB diversity effects

## SB diversity effect on the anoxic to oxic transition for short wait times (1e4)

We showed above that wait times of 1e5 (Figure \@ref(fig:div-temp-meth2)) and 1e4 (Figure \@ref(fig:div-temp-meth3)) produce qualitatively and often quantitatively very similar results to those for wait times of 1e6 (Figure \@ref(fig:div-temp-meth1)). There is, however, a notable difference when there is variation in only SB, or in CB-SB with wait times of 1e4 (Figure \@ref(fig:div-temp-meth3)). Here, moderate increases in diversity have the same effect as with wait times of 1e6. However, larger increases in diversity suddenly cause a large reduction in resilience of the anoxic-oxic transition. Figure \@ref(fig:ss-temp-meth1) gives an example of the stable states when diversity is moderate -- increases in oxygen diffusivity cause strain replacement in the SB to the more tolerant strains, which then delay the shift to the oxic state.  Suprisingly, when there is just a little more diversity in the SB, then the anoxic to oxic transition happens earlier, i.e. there is a large decrease in resilience (Figure \@ref(fig:ss-temp-meth2)).

(ref:ss-temp-meth1) Stable state of the system with **moderate** level of diversity in only SB, and no diversity in the other two groups. All other aspects are as in Figure \@ref(fig:ss-temp-meth3).

```{r ss-temp-meth1, echo = FALSE,  include=TRUE, fig.cap='(ref:ss-temp-meth1)', fig.height=10, fig.width=8, fig.align="center"}
wait_time <- "1e+04"
ss_9s <- readRDS(here::here("data",
                            microxanox_version_path, 
                            "temporal_method",
                            "event_definition_2",
                            paste0("ss_data_9strains_waittime", wait_time, "_event_definition_2.RDS")))
#stab_9s <- readRDS(here::here(\"data", "0_ss_finding", "temporal_method", microxanox_version_path, paste0("stab_data_9strains_waittime", wait_time, "_event_definition_2.RDS"))


#max <- max(ss_3s$PB_var_gmax_s)
#f <- 0.625
#fv <- f*max

only_SB_variable <- ss_9s %>%
  filter(
    abs(CB_var_gmax_s) < 0.0001,
    abs(PB_var_gmax_s) < 0.0001
  )
only_SB_variable <- only_SB_variable[order(only_SB_variable$SB_var_gmax_s), ]

#only_SB_variable$SB_var_gmax_s

fig_state_vs_o2diff_sidebyside_alternative(only_SB_variable[14,])

#fig_state_vs_o2diff(only_SB_variable[13,])
```



(ref:ss-temp-meth2)  Stable state of the system with **higher than moderate** level of diversity in only SB, and no diversity in the other two groups. All other aspects are as in Figure \@ref(fig:ss-temp-meth1).

```{r ss-temp-meth2, echo = FALSE,  include=TRUE, fig.cap='(ref:ss-temp-meth2)', fig.height=10, fig.width=8, fig.align="center"}
fig_state_vs_o2diff_sidebyside_alternative(only_SB_variable[15,])
```

We believe that this phenomenon is caused by an effect that we explain with reference to an analogy of using stepping stones to cross a river. Increasing diversity in our simulations is akin to increasing the width of the river, but keeping the number of stepping stones (strains) constant. So long as the river is not too wide, we can still cross it, but once the width of the river is such that the stepping stones are too far apart, we can't cross the river. And this happens suddenly, the width increases just a tiny bit more, and we suddenly can't cross the river when we could before.

In the simulated microbial ecosystem, the strains are the stepping stones, and it is harder (takes longer) for strains to replace each other when there is greater distance (in trait space) between them. When they can no longer replace each other in the time available, the system behaves as if only the least tolerant, highest growth rate strain (lightest shade of colour) is present. In the case illustrated above, this leads to much lower resilience of the anoxic system state.

This phenomenon is dependent on having static trait values for the strains. If strain traits were dynamic (e.g. due to evolution), then we expect this result to not occur, and rather for the effects of diversity with two or three strains to be the same as those with six or nine.

## SB diversity effect on the oxic to anoxic transition {#SB2}

We show in the main text that diversity in the sulfate-reducing bacteria results in earlier recovery of the anoxic state (Figure 3 c). However, there is only limited evidence of strain replacement on the trajectory of decreasing oxygen diffusivity: the most tolerant strain does increase in density prior to the tipping point but never reaches dominance – it is the strain with lowest tolerance that dominates immediately after the shift to the anoxic state. We here analyze temporal dynamics rather than the final state and show that the more tolerant strains do play a role, but only a transient one. The following dynamics are for oxic starting conditions and with a value of oxygen diffusivity $log_{10}$(oxygen diffusivity) = -4.8. When there is no diversity within functional groups (i.e. two strains that are identical) the system remains oxic (Figure \@ref(fig:sb-div-effect3)). When there is diversity in the two groups of sulfur bacteria the system switches to the anoxic state (Figure \@ref(fig:sb-div-effect4)). (Note that the average trait value [on a log scale] of the two different strains is the same as the trait value of both strains in Figure \@ref(fig:sb-div-effect3).) The more tolerant strain of sulfate-reducing bacteria increases in density in the early phase of the simulation, reduces the concentration of oxygen by inhibiting cyanobacteria, and thus generates an environment in which the less tolerant strain is able to grow. Then, the more tolerant strain is outcompeted by the less tolerant strain as the concentration of the inhibiting substrate (oxygen) declines. 

```{r echo = FALSE}
run_sims_this_section <- TRUE
```

```{r eval = run_sims_this_section, echo = FALSE}
var_expt_levels <- var_expt[,1:6]
no_diversity <- which(rowSums(abs(var_expt_levels))==0)
max_diversty_all <- which(max(rowSums(abs(var_expt_levels))) ==
                            rowSums(abs(var_expt_levels)))
max_only_CB_diversity <- which(max(rowSums(abs(var_expt_levels[,1:2]))) ==
                            rowSums(abs(var_expt_levels[,1:2])) &
                              rowSums(abs(var_expt_levels[,3:6]))==0)
max_only_SBPB_diversity <- which(max(rowSums(abs(var_expt_levels[,3:6]))) ==
                            rowSums(abs(var_expt_levels[,3:6])) &
                              rowSums(abs(var_expt_levels[,1:2]))==0)
```

```{r  eval = run_sims_this_section, echo = FALSE}

num_CB_strains <- 9
num_SB_strains <- 9
num_PB_strains <- 9

sp <- new_strain_parameter(
  n_CB = num_CB_strains,
  n_PB = num_SB_strains,
  n_SB = num_PB_strains,
  values_initial_state = "bush_ssfig3"
)

parameter <- new_runsim_parameter(
  dynamic_model = bushplus_dynamic_model,
  event_definition = event_definition_1,
  event_interval = 100,
  noise_sigma = 0,
  minimum_abundances = rep(1, 3),
  sim_duration = 2000,
  sim_sample_interval = 100,
  strain_parameter = sp,
  log10a_series = c(
    log10(sp$a_O),
    log10(sp$a_O)
  )
)
names(parameter$minimum_abundances) <- c("CB", "PB", "SB")

rm(sp)

parameter$strain_parameter <- var_expt[no_diversity,]$pars[[1]]
parameter$sim_duration <- 100000
parameter$result <- NULL
parameter$log10a_series <- c(-4.8, -4.8)
parameter$strain_parameter$initial_state <- new_initial_state(
  nrow(parameter$strain_parameter$CB),
  nrow(parameter$strain_parameter$PB),
  nrow(parameter$strain_parameter$PB),
  values = "bush_ssfig3"
)
parameter$strain_parameter$initial_state[grep("CB_", names(parameter$strain_parameter$initial_state))] <- 10^10/nrow(parameter$strain_parameter$CB)

sim_res_novar1 <- run_simulation(parameter)

saveRDS(sim_res_novar1, here("data",
                             microxanox_version_path,
                             "temporal_method",
                             "puzzles",
                             "puzzle1_1.RDS"))
```

(ref:sb-div-effect3)  Dynamics of the system **with no diversity** in any functional group, and log10(oxygen diffusivity) = -4.8. The system remains oxic when starting in the oxic state. Although two strains are shown in the legend, they have identical parameters, and so their dynamics are identical and only one is visible in the dynamics.

```{r  sb-div-effect3, echo = FALSE,  include=TRUE, fig.cap='(ref:sb-div-effect3)', fig.height=8, fig.width=8, fig.align="center"}
sim_res_novar1 <- readRDS(here("data",
                             microxanox_version_path,
                             "temporal_method",
                             "puzzles",
                             "puzzle1_1.RDS"))
p <- plot_dynamics(sim_res_novar1)
p[[1]] <- p[[1]] + 
  xlab(NULL) + 
  ylab(expression(log[10](density)~(cells~L^{-1}))) +
  theme_bw()
p[[2]] <- p[[2]] + 
  xlab(NULL) + 
  ylab(expression(log[10](density)~(cells~L^{-1}))) +
  theme_bw()
p[[3]] <- p[[3]] + 
  xlab(NULL) + 
  ylab(expression(log[10](density)~(cells~L^{-1}))) +
  theme_bw()
p[[4]]$data <- p[[4]]$data %>% 
  dplyr::filter(species != "a")
p[[4]]$data$species <- p[[4]]$data$species %>%
  dplyr::recode(SO = "Oxidised sulfur", O = "Oxygen", SR = "Sulfide", P = "Phosphorus")
p[[4]] <- p[[4]] + 
  xlab("time (h)") + 
  ylab(expression(log[10](concentration)~(µM))) + 
  labs(colour = "substrate") +
  theme_bw()
  
p
#ggsave(here("simulations", "expt2", "figures", "switching_novar.pdf"), width = 10)
```


```{r  eval = run_sims_this_section, echo = FALSE}
parameter$strain_parameter <- var_expt[max_only_SBPB_diversity,]$pars[[1]]
#parameter$sim_duration <- 100000
#parameter$result <- NULL
#parameter$log10a_series <- c(-5.2, -5.2)
#parameter$strain_parameter$initial_state <- new_initial_state(
#  nrow(parameter$strain_parameter$CB),
#  nrow(parameter$strain_parameter$PB),
#  nrow(parameter$strain_parameter$PB),
#  values = "bush_ssfig3"
#)
parameter$strain_parameter$initial_state[grep("CB_", names(parameter$strain_parameter$initial_state))] <- 10^10/nrow(parameter$strain_parameter$CB)
sim_res_highvar1 <- run_simulation(parameter)
saveRDS(sim_res_highvar1,
        here("data",
             microxanox_version_path,
             "temporal_method",
             "puzzles",
             "puzzle1_2.RDS"))
```


(ref:sb-div-effect4) Dynamics of the system **with diversity** in the two groups of sulfur bacteria and otherwise exactly the same conditions as in Figure \@ref(fig:sb-div-effect3). The system flips from oxic to anoxic, whereas it did not when there was no diversity.

```{r sb-div-effect4, echo = FALSE,  include=TRUE, fig.cap='(ref:sb-div-effect4)', fig.height=8, fig.width=8, fig.align="center"}
sim_res_highvar1 <- readRDS(here("data",
             microxanox_version_path,
             "temporal_method",
             "puzzles",
             "puzzle1_2.RDS"))

p <- plot_dynamics(sim_res_highvar1)
p[[1]] <- p[[1]] + 
  xlab(NULL) + 
  ylab(expression(log[10](density)~(cells~L^{-1}))) +
  theme_bw()
p[[2]] <- p[[2]] + 
  xlab(NULL) + 
  ylab(expression(log[10](density)~(cells~L^{-1}))) +
  theme_bw()
p[[3]] <- p[[3]] + 
  xlab(NULL) + 
  ylab(expression(log[10](density)~(cells~L^{-1}))) +
  theme_bw()
p[[4]]$data <- p[[4]]$data %>% dplyr::filter(species != "a")
p[[4]]$data$species <- p[[4]]$data$species %>%
  dplyr::recode(SO = "Oxidised sulfur", O = "Oxygen", SR = "Sulfide", P = "Phosphorus")
p[[4]] <- p[[4]] +
  xlab("time (h)") + 
  ylab(expression(log[10](concentration)~(µM))) + 
  labs(colour = "substrate") +
  theme_bw()
  
p

#ggsave(here("simulations", "expt2", "figures", "switching_novar.pdf"), width = 10)
```


# PB diversity effects

## Effect of low to moderate PB diversity

```{r echo = FALSE}
ss_9s <- readRDS(here("data",
                            microxanox_version_path,
                            "temporal_method",
                            "event_definition_2",
                            "ss_data_9strains_waittime1e+06_event_definition_2.RDS"))
only_PB_variable <- ss_9s %>%
  filter(
    abs(CB_var_gmax_s) < 0.0001,
    abs(SB_var_gmax_s) < 0.0001
  )
only_PB_variable <- only_PB_variable[order(only_PB_variable$PB_var_gmax_s), ]

only_PB_stab <- readRDS(here("data",
                                   microxanox_version_path,
                                   "temporal_method",
                                   "event_definition_2",
                                   "processed_data",
                                   "all_stab_data.RDS")) %>%
  filter(
    num_strains == 9,
    waittime == 1e6,
    var_treat == "PB",
    Species == "SB_tot"
  )
only_PB_stab <- only_PB_stab[order(only_PB_stab$PB_var_gmax_s), ]

sim1 <- 2
sim2 <- 10

aO_flip_sim1 <- only_PB_stab[sim1,"hyst_max_log"]
aO_flip_sim2 <- only_PB_stab[sim2,"hyst_max_log"]


```

Increases in the diversity within the PB functional group leads to lower resilience of the anoxic state (Figure \@ref(fig:div-temp-meth1), Panel PB, blue line). For example, a small amount of variation in only the PB functional group, and no variation in other functional groups, results in transition from anoxic to oxic at approximately $log_{10}$(oxygen diffusivity) = `r round(aO_flip_sim1, 2)` (Figure \@ref(fig:pbdiv-only-1)). Whereas a moderate amount of variation in only the PB functional group, and no variation in other functional groups, results in transition from anoxic to oxic at approximately $log_{10}$(oxygen diffusivity) = `r round(aO_flip_sim2, 2)` (Figure \@ref(fig:pbdiv-only-2)).


(ref:pbdiv-only-1) Right-hand panels show the transition from anoxic to oxic state when oxygen diffusivity decreases. The results here are for a **low amount of variation in only the PB functional group**. The transition from anoxic to oxic occurs at approximately $log_{10}$(oxygen diffusivity) = `r round(aO_flip_sim1, 2)`. All other aspects are the same as in Figure \@ref(fig:ss-temp-meth3).

```{r pbdiv-only-1, echo = FALSE,  include=TRUE, fig.cap='(ref:pbdiv-only-1)', fig.height=10, fig.width=8, fig.align="center"}
fig_state_vs_o2diff_sidebyside_alternative(only_PB_variable[sim1,])
```

(ref:pbdiv-only-2) Right-hand panels show the transition from anoxic to oxic state when oxygen diffusivity decreases. The results here are for a **moderate amount of variation in only the PB functional group**. The transition from anoxic to oxic occurs at approximately $log_{10}$(oxygen diffusivity) = `r round(aO_flip_sim2, 2)`. All other aspects are the same as in Figure \@ref(fig:ss-temp-meth3).

```{r pbdiv-only-2, echo = FALSE,  include=TRUE, fig.cap='(ref:pbdiv-only-2)', fig.height=10, fig.width=8, fig.align="center"}
fig_state_vs_o2diff_sidebyside_alternative(only_PB_variable[sim2,])
```


This result is different to the effect of diversity within the CB and SB functional groups, which both cause increases in resilience of the state that they dominate (Figure \@ref(fig:div-temp-meth1)). Why does diversity in the PB functional group reduce the resilience of the state it dominates?

The counter-intuitive effect of diversity in the phototrophic sulfur bacteria (PB) possibly results from their positive effects on conditions and processes that favour the oxic state. The phototrophic sulfur bacteria oxidise reduced sulfur to oxidised sulfur. Since reduced sulfur inhibits the cyanobacteria (CB), PB can have a positive effect on CB via their consumption of reduced sulfur. Hence, high growth rates of PB will have a positive effect on CB, and therefore favour the oxic state. We would expect to see that higher diversity of PB is associated with lower concentrations of reduced sulfur and higher concentrations of oxidised sulfur before the transition.

Another pathway by which a higher growth rate of PB can favour the oxic state is by indirectly increasing the concentration of oxygen, thus inhibiting themselves and sulfate-reducing bacteria. This effect occurs because PB consume reduced sulfur, such that less reduced sulfur is available for abiotic oxidation. Consequently, a larger amount of oxygen remains in the environment and inhibits sulfur bacteria. We would expect to see that higher diversity of PB is associated with higher concentrations of oxygen before the transition. 

Among these two pathways, effects of phototrophic sulfur bacteria via abiotic oxidation of sulfide are probably of subordinate importance in natural environments because abiotic oxidation rates of sulfide are very small compared to biotic oxidation rates ([Luther et al 2011]( https://doi.org/10.3389/fmicb.2011.00062)). 

There is, however, at least one pathway by which PB can favour the anoxic state: PB produce oxidised sulfur which is consumed by SB, and SB produce reduced sulfur, which inhibits CB.

Looking at the strain replacement that occurs across the oxygen diffusivity gradient, we see that as oxygen diffusivity increases, the more oxygen tolerant strains of PB replace the less tolerant ones because they have higher realised growth rates under the prevailing environmental conditions. Greater PB diversity leads to the availability of strains with higher oxygen tolerance and thus even higher realised growth rates, which via the pathways described above, can favour the oxic state.

Substrate concentrations close to the anoxic-oxic transition in simulations with low and high PB diversity, respectively, corroborate the expected effects (Figure \@ref(fig:pbdiv-only-3)). As the transition is approached, higher PB diversity leads to higher concentration of oxygen, higher concentration of oxidised sulfur (SO), and lower concentration of reduced sulfur (SR), than a simulation with lower PB diversity. All these differences in substrate concentrations are in the direction of favouring the oxic state, resulting in the earlier transition to the oxic state.

(ref:pbdiv-only-3) Concentrations of three substrates (O: oxygen; SO: oxidised sulfur; SR: reduced sulfur) in a low PB diversity simulation (x-axis) versus substrate concentration in a moderate PB diversity simulation (y-axis). (The same levels of PB diversity as in Figures \@ref(fig:pbdiv-only-1) and \@ref(fig:pbdiv-only-2), respectively). The figures show only a small range of oxygen diffusivity on the trajectory of increasing oxygen diffusivity. The black line is the 1:1 line.


```{r pbdiv-only-3, echo = FALSE,  include=TRUE, fig.cap='(ref:pbdiv-only-3)', fig.height=5, fig.width=8, fig.align="center"}


stab_9s <- readRDS(here("data", 
                        microxanox_version_path,
                        "temporal_method", 
                        "event_definition_2",
                        "ss_data_9strains_waittime1e+06_event_definition_2.RDS"))



ss1 <- only_PB_variable[sim1,]
ss2 <- only_PB_variable[sim2,]

dd1 <- ss1$ssfind_result[[1]] %>%
  pivot_longer(names_to = "State_variable", values_to = "Quantity_sim1", 2:32)
dd2 <- ss2$ssfind_result[[1]] %>%
  pivot_longer(names_to = "State_variable", values_to = "Quantity_sim2", 2:32)
dd <- full_join(dd1, dd2)
rm(ss1, ss2, dd1, dd2)

labs <- seq(-3.5, -3.3, by=0.05)

dd$State_variable <- dd$State_variable %>%
  dplyr::recode(SO = "Oxidised sulfur", O = "Oxygen", SR = "Sulfide", P = "Phosphorus")


dd %>%
  filter(State_variable %in% c("Oxygen", "Sulfide", "Oxidised sulfur"),
         direction == "up",
         a_O < as.numeric(aO_flip_sim2),
         a_O > -3.5) %>%
  ggplot(aes(x = log10(Quantity_sim1),
                 y = log10(Quantity_sim2),
                 col = a_O )) +
  geom_point(lwd = 3) +
  geom_path() +
  scale_color_gradient(low = "yellow", high = "red", na.value = NA,
                       breaks = labs,
                       labels = labs) +
  #coord_trans(x = "log10", y = "log10") +
  facet_wrap( ~ State_variable, scales = "free") +
  geom_abline(intercept = 0, slope = 1) +
  xlab(expression(atop(log[10](concentration)~(µM), Low~PB~diversity))) +
  ylab(expression(atop(log[10](concentration)~(µM),Moderate~PB~diversity))) +
  labs(col=expression(paste("log"[10],"(oxygen diffusivity) (h"^{-1}, ")"))) +
  theme_bw() +
  theme(legend.position="top")
```


## Effect of medium to high PB diversity


When there was medium to high variation in only the phototrophic sulfur bacteria (PB), the strain dynamics deviated from the standard pattern: at low oxygen diffusivity, all three functional groups coexisted. Specifically, on the trajectory of decreasing oxygen diffusivity, cyanobacteria reappeared at low oxygen diffusivity after their collapse (medium variation in PB, Figure \@ref(fig:PB-stabstat1)) or never fully collapsed (high variation in PB, Figure \@ref(fig:PB-stabstat2)). On the trajectory of increasing oxygen diffusivity, cyanobacteria coexisted with the two groups of sulfur bacteria at low oxygen diffusivity, collapsed at medium levels of oxygen diffusivity, and then increased to sole dominance at the tipping point to the oxic state (Figure \@ref(fig:PB-stabstat1), Figure \@ref(fig:PB-stabstat2)). 

The same pattern occurred for medium to high simultaneous variation in phototrophic sulfur bacteria and cyanobacteria, albeit only on the trajectory of increasing oxygen diffusivity (Figure \@ref(fig:CBPB-stabstat3)). On the trajectory of decreasing oxygen diffusivity, cyanobacteria dominated over the entire gradient of oxygen diffusivity, such that effects of PB diversity did not play out on this trajectory. 

The reason for coexistence of all three groups at low oxygen diffusivity is likely that the PB strain that dominated at low oxygen diffusivity (the least tolerant strain) lowered the sulfide concentrations to such low levels that cyanobacteria were able to grow. The higher the trait variation in PB, the higher the maximum growth rate of the least tolerant PB strain, and the higher its consumption of sulfide. Consequently, moderate variation in PB lowered the sulfide concentration such that the tipping point from anoxic to oxic occurred earlier (as discussed in the previous section), and medium to high variation in PB lowered the sulfide concentration even further such that cyanobacteria were able to grow and coexist with the two groups of sulfur bacteria at low oxygen diffusivity. 

At intermediate levels of oxygen diffusivities, the cyanobacteria were excluded (Figure \@ref(fig:PB-stabstat1)). We speculate that this might be associated with the increase in sulfide prior to the collapse of the phototrophic sulfur bacteria. This increase in sulfide occurred simultaneously with turnover to PB strains with lower maximum growth rate and with a decline in PB density; these two processes likely resulted in less consumption of sulfide. 


```{r}
wait_time <- "1e+06"
ss_9s <- readRDS(here::here("data",
                            microxanox_version_path, 
                            "temporal_method",
                            "event_definition_2",
                            paste0("ss_data_9strains_waittime", wait_time, "_event_definition_2.RDS")))



max_CB_var_gmax <- max(ss_9s$CB_var_gmax_s)
max_SB_var_gmax <- max(ss_9s$SB_var_gmax_s)
max_PB_var_gmax <- max(ss_9s$PB_var_gmax_s)

ss_9s <- ss_9s %>%
  mutate(standvar_CB_gmax = CB_var_gmax_s / max_CB_var_gmax,
         standvar_SB_gmax = SB_var_gmax_s / max_SB_var_gmax,
         standvar_PB_gmax = PB_var_gmax_s / max_PB_var_gmax)
```


(ref:PB-stabstat1) Coexistence of all three groups at low oxygen diffusivities when only PB contains strain variation, and the amount of strain variation is medium to high. This example is when the standardised amount of trait variation is 0.74 (0 being no variation, and 1 being the maximum amount in any simulation.)

```{r PB-stabstat1, echo = FALSE,  include=TRUE, fig.cap='(ref:PB-stabstat1)', fig.height=10, fig.width=8, fig.align="center"}
this_dat <- ss_9s %>%
  filter(
    CB_var_gmax_s == 0,
    SB_var_gmax_s == 0
  )
this_dat <- this_dat[order(this_dat$PB_var_gmax_s), ]
which_one <- 15
#print(paste0("Standardised level of variation = ",
#round(this_dat$standvar_PB_gmax[this_one],2)))
fig_state_vs_o2diff_sidebyside_alternative(this_dat[which_one,])
```



(ref:PB-stabstat2) High trait variation in only the PB group results in coexistence of all groups at low oxygen diffusivity levels (as in Figure \@ref(fig:PB-stabstat1)) and also no collapse in CB abundance in the decreasing oxygen diffusivity trajectory. This example is when the standardised amount of trait variation is 0.95 (0 being no variation, and 1 being the maximum amount in any simulation.)

```{r PB-stabstat2, echo = FALSE,  include=TRUE, fig.cap='(ref:PB-stabstat2)', fig.height=10, fig.width=8, fig.align="center"}
which_one <- 19
#print(paste0("Standardised level of variation = ",
#round(this_dat$standvar_PB_gmax[this_one],2)))
fig_state_vs_o2diff_sidebyside_alternative(this_dat[which_one,])
```



(ref:CBPB-stabstat3) Stable states when there is high trait variation in both CB and PB groups. This example is when the standardised amount of trait variation is 0.95 (0 being no variation, and 1 being the maximum amount in any simulation.)


```{r CBPB-stabstat3, echo = FALSE,  include=TRUE, fig.cap='(ref:CBPB-stabstat3)', fig.height=10, fig.width=8, fig.align="center"}
this_dat <- ss_9s %>%
  filter(
    SB_var_gmax_s == 0,
    abs(standvar_CB_gmax - standvar_PB_gmax) < 0.00001
  )
this_dat <- this_dat[order(this_dat$CB_var_gmax_s), ]
which_one <- 19
#print(paste0("Standardised level of variation = ",
#round(this_dat$standvar_PB_gmax[this_one],2)))
fig_state_vs_o2diff_sidebyside_alternative_leftfliponly(this_dat[which_one,])
```





# Underlying mechanisms

In the main text we show that trait variation (functional diversity) can affect resilience. In this supplementary section, we discuss the mechanism(s) underlying the observed diversity effects. We first explain the two major classes of mechanisms known to drive diversity effects, and then use additional simulations to elucidate the mechanisms operating in our study.

## Complementarity and sampling effect

Research about the relationship between biodiversity and ecosystem functioning has identified two major mechanisms by which diversity affects ecosystem functioning: the selection/sampling effect and the complementarity effect (Loreau & Hector 2001). The sampling effect occurs when more diverse communities have higher functioning because they have a higher probability of containing a high performing species, and when only one (or few) of the species in a community determine the community and ecosystem level properties, such as biomass production or stability. In this case the observed effect of diversity is driven by only one (or few) species (Cardinale et al. 2011).  In contrast, the complementarity effect occurs when higher ecosystem functioning in more diverse communities is due to niche differences or positive interactions among species, i.e. the diversity effect is driven by combinations of the species in a community (Cardinale et al. 2011). Both effects are considered valid and important mechanisms by which diversity affects ecosystem functioning (Cardinale et al. 2006, 2011, 2012; Fargione et al. 2007).

Similarly, the complementarity effect and the sampling effect are important mechanisms by which diversity affects stability (Loreau & Mazancourt 2013; Loreau et al. 2021) and resilience (Steiner et al. 2006; Hughes & Stachowicz 2011). For example, diversity can stabilize aggregate ecosystem properties (e.g. community biomass) when species differ in their responses to environmental factors and therefore respond asynchronously to environmental fluctuations (temporal complementarity effect; Loreau & Mazancourt 2013). However, the sampling effect might be of particular importance when considering ecosystem responses to directional environmental change or pulse disturbances rather than to environmental fluctuations (Steiner et al. 2006): more diverse systems might absorb a greater amount of change by having a higher chance of containing species with high environmental tolerance. In this scenario, the positive diversity effect on resilience is driven by one (or few) species. From a conservation perspective, however, one would nevertheless strive to maintain diversity and not only the most tolerant species because of the uncertainties associated with environmental change: we do not know with sufficient certainty which stressors will change in a given ecosystem and which traits will therefore be relevant in the future.




```{r samp-effect1, echo = FALSE, eval= FALSE, include=TRUE, fig.cap='(ref:samp-effect1)', fig.height=5, fig.width=6, fig.align="center"}

## no longer included in supplement.

## (ref:samp-effect1) A simple phenomenological model of the effect of number of strains in an ecosystem on the resilience of that ecosystem, when the resilience of a ecosystem is linearly proportional to the maximum trait value of the strains present in that ecosystem. Each data point is a different ecosystem, with the strains in each ecosystem drawn at random from a pool of 20 strains. Each of the 20 strains was assigned a trait value drawn at random from a uniform distribution with minimum 0 and maximum 1. The blue line is the fit of a generalised additive model and the grey ribbon shows the 1 standard error confidence range.


if(!("tidyverse" %in% .packages()))
  library(tidyverse)
set.seed(12345)
number_of_species <- 20
species_trait <- tibble(species_name = paste0("Species-", letters[1:20]),
                        trait_value = runif(number_of_species))

richness_levels <- seq(1, 10, 1)
compositions_per_richness_level <- 10
expt <- crossing(richness = richness_levels,
                 composition = paste0("Comp-", 1:compositions_per_richness_level)) %>%
  mutate(composition = paste0(composition, "-", richness),
         resilience_max_trait = NA)

for(i in 1:nrow(expt))
  expt$resilience_max_trait[i] = max(sample(species_trait$trait_value, expt$richness[i]))

ggplot(expt, aes(x = richness, y = resilience_max_trait)) +
  geom_point() +
  xlab("Number of strains") +
  ylab("Resilience if only the tolerance of\nthe most tolerant strain determines resilience.") +
  geom_smooth(method = "gam")


```

## Investigating the mechanisms operating

We investigated if the diversity effects observed in our study were in line with the sampling effect, that is, driven entirely by the presence of the most tolerant strain or if other strains played a role as well. To this end, we ran additional simulations where functional groups contained only the strain with greatest tolerance (which due to the tradeoff also had the lowest maximum growth rate). That is, rather than increasing diversity we only increased tolerance of a single strain. Specifically, for each of the 20 levels of trait variation, we made simulations where each functional group had only one strain, and its trait values corresponded with those of the most tolerant strain for the given level of trait variation. We also manipulated in which of the functional groups tolerance increased. For example, in a simulation to be compared with variation in only the SB group (and no variation in the CB or PB groups), the SB group contained the strain with highest tolerance, whereas CB and PB each contained the strain with mean tolerance. We then compared the tipping point positions resulting from this tolerance manipulation with the tipping point positions resulting from the diversity manipulation presented in the main text. In addition, we compared the strain dynamics across the oxygen diffusivity gradient of the tolerance manipulation and the diversity manipulation. If the observed diversity effects are caused by the sampling effect, we expect there to be no difference between the tolerance manipulation and the diversity manipulation. 

We found no difference between the diversity and the tolerance manipulation when trait values were low (Figure \@ref(fig:max-tol1)). When trait values were medium to high, however, we did find differences between the diversity and tolerance manipulation for most combinations except for diversity/tolerance in only CB and in SB-PB (Figure \@ref(fig:max-tol1)). Specifically, the diversity and tolerance manipulation differed either in tipping point positions (CB-SB and CB-SB-PB) or in whether coexistence patterns were atypical (SB, PB, CB-PB, CB-SB). We describe and discuss these differences in detail in the sections below. Briefly, in these scenarios less tolerant strains did play a role because of their high maximum growth rates. For example, absence of SB strains with high growth rate led to reduced production of sulfide and therefore to coexistence with cyanobacteria at low oxygen diffusivity or to increased resilience of the oxic state. In contrast, absence of PB strains with high growth rate led to reduced consumption of sulfide and therefore precluded coexistence with cyanobacteria at low oxygen diffusivity (as observed when there was high variation in PB). 

Taken together, at low amounts of trait variation the diversity effects were driven entirely by the most tolerant strain, in line with a sampling effect. At medium to high trait variation, however, in most combinations both tolerant and less tolerant strains contributed to the diversity effects on tipping point positions or coexistence patterns. The reason was that in the two groups of sulfur bacteria the trade-off between tolerance and growth rate had ecosystem-level effects, i.e. effects on the sulfide concentration and therefore on the cyanobacteria. 


(ref:max-tol1) Effects of strain variation on position of tipping points compared between simulations with all (nine) strains present and simulations where only the most tolerant strain is present. When tipping points are shown with dashed lines then patterns of coexistence were atypical of the classical pattern shown in Figure 1a of the main text (e.g. when cyanobacteria coexisted with sulfur bacteria at low oxygen diffusivities).

```{r max-tol1, echo = FALSE,  include=TRUE, fig.cap='(ref:max-tol1)', fig.height=8, fig.width=8, fig.align="center"}
#all_stab <- readRDS(here("data", "0_ss_finding", "temporal_method", "processed_data", "stab_data_temporal_method_ELErevision1.2.RDS"))

## read the data
all_stab <- readRDS(here("data",
                         microxanox_version_path,
                         "temporal_method", 
                         "event_definition_2",
                         "processed_data",
                         "all_stab_data.RDS"))

which_strains <- c(-1, 9)
which_wait_time <- 1e6
figure_title <- " "


plot_here11 <- all_stab %>%
  mutate(stab_states_type = "standard",
         stab_states_type = case_when(
          (num_strains == 9 & var_treat == "CB" & stand_var > 1.00) ~ "non-standard",
          (num_strains == 9 & var_treat == "SB" & stand_var > 1.00) ~ "non-standard",
          (num_strains == 9 & var_treat == "PB" & stand_var > 0.65) ~ "non-standard",
          (num_strains == 9 & var_treat == "CB-SB" & stand_var > 1.00) ~ "non-standard",
          (num_strains == 9 & var_treat == "CB-PB" & stand_var > 0.65) ~ "non-standard",
          (num_strains == 9 & var_treat == "SB-PB" & stand_var > 1.00) ~ "non-standard",
          (num_strains == 9 & var_treat == "CB-SB-PB" & stand_var > 1.00) ~ "non-standard",
          
          (num_strains == -1 & var_treat == "CB" & stand_var > 1.00) ~ "non-standard",
          (num_strains == -1 & var_treat == "SB" & stand_var > 0.60) ~ "non-standard",
          (num_strains == -1 & var_treat == "PB" & stand_var > 1.00) ~ "non-standard",
          (num_strains == -1 & var_treat == "CB-SB" & stand_var > 0.60) ~ "non-standard",
          (num_strains == -1 & var_treat == "CB-PB" & stand_var > 1.00) ~ "non-standard",
          (num_strains == -1 & var_treat == "SB-PB" & stand_var > 1.00) ~ "non-standard",
          (num_strains == -1 & var_treat == "CB-SB-PB" & stand_var > 1.00) ~ "non-standard"
         ),
         stab_states_type = ifelse(is.na(stab_states_type), "standard", "non-standard")
  ) %>%
  mutate(hyst_min_log = ifelse(var_treat == "CB-SB" &
                                 stab_states_type == "non-standard" &
                                 Species == "SB_tot" &
                                 stand_var > 0.85,
                               -8,
                               hyst_min_log
                                 ),
         hyst_max_log = ifelse(var_treat == "CB-SB" &
                                 stab_states_type == "non-standard" &
                                 Species == "SB_tot" &
                                 stand_var > 0.85,
                               -8,
                               hyst_max_log
                                 ))
   
plot_1 <- fig_div_vs_o2diff_maxtolcomp_stable_state_type(plot_here11, which_strains,
                              which_wait_time, figure_title)
plot_1 +
  xlab("Standardised amount of trait variation used to create nine strains\nof which either all nine are retained, or only the most tolerant is retained.") +
  guides(col=guide_legend(title="Which strains are retained", title.position = "top"))




```


```{r echo = FALSE}
nine_strains <- readRDS(here("data",
                         microxanox_version_path,
                         "temporal_method", 
                         "event_definition_2",
                          "ss_data_9strains_waittime1e+06_event_definition_2.RDS"))
maxtolonly_strains <- readRDS(here("data",
                         microxanox_version_path,
                         "temporal_method", 
                         "event_definition_2",
                          "ss_data_-1strains_waittime1e+06_event_definition_2.RDS"))

#nine_strains$parameter[[1]]$strain_parameter
#maxtolonly_strains[134,]$parameter[[1]]$strain_parameter
```

## Increasing tolerance vs. increasing diversity in only the SB

When SB contained multiple strains varying in growth rate (and tolerance) there were two stable states at intermediate levels of oxygen diffusivity, only one stable state at low oxygen diffusivity with only sulfur bacteria, and only one stable state at high oxygen diffusivity with high cyanobacterial density (Figure \@ref(fig:max-tol2)), as was previously reported. When, however, only one SB strain with very high tolerance and low maximum growth rate was present, the stable states change: at low oxygen diffusivity  all three functional groups coexisted (Figure \@ref(fig:max-tol3)) – this never occurred when there were multiple SB strains varying in tolerance.

Apparently, at low oxygen diffusivity levels, the absence of SB strains with high growth rate precludes the formation of a stable state that does not include cyanobacteria. The most tolerant SB strain has a low maximum growth rate, which leads to less production of sulfide. Presence of only the most tolerant SB strain therefore results in lower sulfide concentrations (Figure \@ref(fig:max-tol2) vs Figure \@ref(fig:max-tol3)), allowing growth of cyanobacteria at low oxygen diffusivity.

However, cyanobacteria are excluded at intermediate oxygen diffusivity. We speculate that this might be associated with the increase in sulfide at intermediate oxygen diffusivity, which might be due to the decline and collapse of PB and therefore reduced consumption of sulfide.

(ref:max-tol2) Stable states found when there is variation among nine strains of SB and the amount of variation is high. All other aspects are the same as in Figure \@ref(fig:ss-temp-meth3). 

```{r max-tol2, echo = FALSE,  include=TRUE, fig.cap='(ref:max-tol2)', fig.height=10, fig.width=8, fig.align="center"}
## max SB_var_gmax_s is 0.09473684
nine_strain_one <- filter(nine_strains,
       CB_var_gmax_s == 0,
       abs(SB_var_gmax_s - 0.09473684) < 0.001,
       PB_var_gmax_s == 0)
ss_result <- nine_strain_one
fig_state_vs_o2diff_sidebyside_alternative_leftfliponly(nine_strain_one)
```

(ref:max-tol3) Stable states found when SB contained only the most tolerant strain, with trait values corresponding with the diversity level of the previous figure. All other aspects are the same as in Figure \@ref(fig:ss-temp-meth3). 

```{r max-tol3, echo = FALSE,  include=TRUE, fig.cap='(ref:max-tol3)', fig.height=10, fig.width=8, fig.align="center"}
max_tol_one <- filter(maxtolonly_strains,
       CB_var_gmax_s == 0,
       abs(SB_var_gmax_s - 0.09473684) < 0.001,
       PB_var_gmax_s == 0)
fig_state_vs_o2diff_sidebyside_alternative(max_tol_one)
```

In the ecosystem shown in Figure \@ref(fig:max-tol3), where the SB contained one strain with very high tolerance, there was large variation in state variables along the oxygen diffusivity gradient at values less than about -4.5. To examine what might be causing this, we simulated the ecosystem at a constant oxygen diffusivity of $log_{10}$ oxygen diffusivity of -6 (Figure \@ref(fig:SB-maxtol-dyn1)). We see that the dynamics are qualitatively different from other cases in which there is a stable state reached in such a simulation (Figure \@ref(fig:sb-div-effect4)). It appears that the high tolerance SB strain is able to initiate a transition to an anoxic state, but is not able to create a stable anoxic state. This dynamical situation was not present in the ecosystem in which high tolerance and low tolerance strains were present.

```{r echo = FALSE}
#nine_strain_one$pars[[1]]
#max_tol_one$pars[[1]]
#fig_state_vs_o2diff_sidebyside_alternative(nine_strain_one)
#fig_state_vs_o2diff_sidebyside_alternative(max_tol_one)
parameter <- max_tol_one$parameter[[1]]
parameter$strain_parameter <- max_tol_one$pars[[1]]
#parameter$strain_parameter <- ss_result$pars[[1]]
total_initial_abundances <- 10^5
parameter <- set_temporal_ssfind_initial_state(parameter,
                                               total_initial_abundances,
                                               total_initial_abundances,
                                               total_initial_abundances)
parameter$log10a_series <- c(-6, -6)
wait_time <-  1e5
parameter$sim_duration <- wait_time * length(parameter$log10a_series)
parameter$sim_sample_interval <- 10
#parameter$sim_duration
test1 <- run_simulation(parameter)
```

(ref:SB-maxtol-dyn1) Temporal dynamics of an ecosystem with a single SB strain with high tolerance, the same one as in Figure \@ref(fig:max-tol3) simulated with constant $log_{10}$ oxygen diffusivity of -6.

```{r SB-maxtol-dyn1, echo = FALSE,  include=TRUE, fig.cap='(ref:SB-maxtol-dyn1)', fig.height=10, fig.width=8, fig.align="center"}
p <- plot_dynamics(test1)
p[[1]] <- p[[1]] + 
  xlab(NULL) + 
  ylab(expression(log[10](density)~(cells~L^{-1}))) +
  theme_bw()
p[[2]] <- p[[2]] + 
  xlab(NULL) + 
  ylab(expression(log[10](density)~(cells~L^{-1}))) +
  theme_bw()
p[[3]] <- p[[3]] + 
  xlab(NULL) + 
  ylab(expression(log[10](density)~(cells~L^{-1}))) +
  theme_bw()
p[[4]]$data <- p[[4]]$data %>% 
  dplyr::filter(species != "a")
p[[4]]$data$species <- p[[4]]$data$species %>%
  dplyr::recode(SO = "Oxidised sulfur", O = "Oxygen", SR = "Sulfide", P = "Phosphorus")
p[[4]] <- p[[4]] + 
  xlab("time (h)") + 
  ylab(expression(log[10](concentration)~(µM))) + 
  labs(colour = "substrate") +
  theme_bw()
  
p
```




## Increasing tolerance vs. increasing diversity in CB-SB

In the CB-SB case, the presence of only the most tolerant strain led to coexistence of all groups at low levels of oxygen diffusivity when trait values were moderately high (Figure \@ref(fig:max-tol55)) or to absence of alternate stable states and sole dominance of cyanobacteria when trait values were high (Figure \@ref(fig:max-tol5)). As discussed in the previous section, the low maximum growth rate of the most tolerant SB strain led to less production of sulfide such that SB could not exclude CB. 

(ref:max-tol55)  Stable states found when CB and SB each contain only the single most tolerant strain with moderate tolerance. All other aspects are the same as in Figure \@ref(fig:ss-temp-meth3). 

```{r max-tol55, echo = FALSE,  include=TRUE, fig.cap='(ref:max-tol55)', fig.height=10, fig.width=8, fig.align="center"}
#max_tol_one_CBSB <- filter(maxtolonly_strains,
#       abs(CB_var_gmax_s - 0.026592798) < 0.001,
#       abs(SB_var_gmax_s - 0.079778395) < 0.001,
#       PB_var_gmax_s == 0)
max_tol_one_CBSB <- filter(maxtolonly_strains,
       abs(CB_var_gmax_s - 0.023268699) < 0.001,
       abs(SB_var_gmax_s - 0.069806096) < 0.001,
       PB_var_gmax_s == 0)


fig_state_vs_o2diff_sidebyside_alternative(max_tol_one_CBSB)
```



(ref:max-tol5)  Stable states found when CB and SB each contain only the single most tolerant strain with high tolerance. All other aspects are the same as in Figure \@ref(fig:ss-temp-meth3). 

```{r max-tol5, echo = FALSE,  include=TRUE, fig.cap='(ref:max-tol5)', fig.height=10, fig.width=8, fig.align="center"}
max_tol_one_CBSB <- filter(maxtolonly_strains,
       abs(CB_var_gmax_s - 0.026592798) < 0.001,
       abs(SB_var_gmax_s - 0.079778395) < 0.001,
       PB_var_gmax_s == 0)
max_tol_one_CBSB <- filter(maxtolonly_strains,
       abs(CB_var_gmax_s - 0.028254848) < 0.001,
       abs(SB_var_gmax_s - 0.084764545) < 0.001,
       PB_var_gmax_s == 0)


fig_state_vs_o2diff_sidebyside_alternative_nofliparrows(max_tol_one_CBSB)
```



## Increasing tolerance vs. increasing diversity in PB and in PB-CB

Whereas medium to high diversity in only PB led to coexistence of all three functional groups at low oxygen diffusivity levels (Figure \@ref(fig:PB-stabstat4)), presence of only the most tolerant PB strain led to the standard-pattern of either dominance of sulfur bacteria or of cyanobacteria (Figure \@ref(fig:PB-stabstat5)). The most tolerant PB strain has a low maximum growth rate, leading to less consumption of sulfide and therefore to higher sulfide concentrations compared to a system that contains both slow- and fast-growing strains. Due to the higher sulfide concentrations in the system with only the most tolerant PB strain, cyanobacteria were not able to grow at low oxygen diffusivity levels. 

For the same reasons we found a similar pattern on the trajectory of increasing oxygen diffusivity for medium to high simultaneous trait variation in PB and CB but not for the corresponding tolerance manipulation (data not shown).  

```{r}
wait_time <- "1e+06"
ss_9s <- readRDS(here::here("data",
                            microxanox_version_path, 
                            "temporal_method",
                            "event_definition_2",
                            paste0("ss_data_9strains_waittime", wait_time, "_event_definition_2.RDS")))



max_CB_var_gmax <- max(ss_9s$CB_var_gmax_s)
max_SB_var_gmax <- max(ss_9s$SB_var_gmax_s)
max_PB_var_gmax <- max(ss_9s$PB_var_gmax_s)

ss_9s <- ss_9s %>%
  mutate(standvar_CB_gmax = CB_var_gmax_s / max_CB_var_gmax,
         standvar_SB_gmax = SB_var_gmax_s / max_SB_var_gmax,
         standvar_PB_gmax = PB_var_gmax_s / max_PB_var_gmax)
```


(ref:PB-stabstat4) Coexistence of all three groups at low oxygen diffusivities when only PB contains strain variation, and the amount of strain variation is medium to high. This example is when the standardised amount of trait variation is 0.89 (0 being no variation, and 1 being the maximum amount in any simulation.)

```{r PB-stabstat4, echo = FALSE,  include=TRUE, fig.cap='(ref:PB-stabstat4)', fig.height=10, fig.width=8, fig.align="center"}
this_dat <- ss_9s %>%
  filter(
    CB_var_gmax_s == 0,
    SB_var_gmax_s == 0
  )
this_dat <- this_dat[order(this_dat$PB_var_gmax_s), ]
which_one <- 18
#print(paste0("Standardised level of variation = ",
#round(this_dat$standvar_PB_gmax[this_one],2)))
fig_state_vs_o2diff_sidebyside_alternative(this_dat[which_one,])
```


```{r}
wait_time <- "1e+06"
ss_tol <- readRDS(here::here("data",
                            microxanox_version_path, 
                            "temporal_method",
                            "event_definition_2",
                            paste0("ss_data_-1strains_waittime", wait_time,
                                   "_event_definition_2.RDS")))



max_CB_var_gmax <- max(ss_tol$CB_var_gmax_s)
max_SB_var_gmax <- max(ss_tol$SB_var_gmax_s)
max_PB_var_gmax <- max(ss_tol$PB_var_gmax_s)

ss_tol <- ss_tol %>%
  mutate(standvar_CB_gmax = CB_var_gmax_s / max_CB_var_gmax,
         standvar_SB_gmax = SB_var_gmax_s / max_SB_var_gmax,
         standvar_PB_gmax = PB_var_gmax_s / max_PB_var_gmax)
```


(ref:PB-stabstat5) Dynamics when PB contain only the most tolerant strain, with this being the same tolerance level as the most tolerant strain in Figure \@ref(fig:PB-stabstat4) immediately above.

```{r PB-stabstat5, echo = FALSE,  include=TRUE, fig.cap='(ref:PB-stabstat5)', fig.height=10, fig.width=8, fig.align="center"}
this_dat <- ss_tol %>%
  filter(
    CB_var_gmax_s == 0,
    SB_var_gmax_s == 0
  )
this_dat <- this_dat[order(this_dat$PB_var_gmax_s), ]
which_one <- 18
#print(paste0("Standardised level of variation = ",
#round(this_dat$standvar_PB_gmax[this_one],2)))
fig_state_vs_o2diff_sidebyside_alternative(this_dat[which_one,])
```




## Increasing tolerance vs. increasing diversity in CB-SB-PB

At high levels of diversity / high levels of tolerance, resilience of the oxic state is greater in the tolerance than diversity manipulation. We speculate that this is because the most tolerant SB strain has reduced ability to suppress the cyanobacteria, as we have also observed in the scenarios of tolerance in only SB and in SB-CB (described above). 

(ref:max-tol6) Stable states found with nine strains of CB-SB-PB and high diversity. All other aspects are the same as in Figure \@ref(fig:ss-temp-meth3). 

```{r max-tol6, echo = FALSE,  include=TRUE, fig.cap='(ref:max-tol6)', fig.height=10, fig.width=8, fig.align="center"}
filter(nine_strains,
       abs(CB_var_gmax_s - 0.026592798) < 0.001,
       abs(SB_var_gmax_s - 0.079778395) < 0.001,
       abs(PB_var_gmax_s - 0.079778395) < 0.001) %>%
  fig_state_vs_o2diff_sidebyside_alternative()

```

(ref:max-tol7) Stable states found with the single most tolerant CB strain, SB strain, and PB strain, with trait values corresponding with the diversity level in the previous figure. All other aspects are the same as in Figure \@ref(fig:ss-temp-meth3). 

```{r max-tol7, echo = FALSE,  include=TRUE, fig.cap='(ref:max-tol7)', fig.height=10, fig.width=8, fig.align="center"}
filter(maxtolonly_strains,
       abs(CB_var_gmax_s - 0.026592798) < 0.001,
       abs(SB_var_gmax_s - 0.079778395) < 0.001,
       abs(PB_var_gmax_s - 0.079778395) < 0.001) %>%
  fig_state_vs_o2diff_sidebyside_alternative()
```




# (Non-)additivity

We investigated if the effects of trait variation of individual functional groups were additive. We compared the tipping point positions as determined by our simulations with the expected position of the respective tipping point if effects of diversity in different functional groups were additive. The expected tipping point was calculated by adding the effect sizes (i.e. the difference in tipping point position between the no-diversity scenario and a given amount of trait variation) on a linear scale. 


```{r fig-5-log, echo = FALSE,  include=TRUE, fig.cap='Previously figure 5 Main Manuscript log', fig.height=8, fig.width=8, fig.align="center"}
## Figure 5 ----
## Assessing additivity of diversity effect sizes
## By Uriah
## Owen switched to SB_tot

all_stab_results_small <- readRDS(here("data",
                         microxanox_version_path,
                         "temporal_method", 
                         "event_definition_2",
                         "processed_data", "all_stab_data.RDS"))

wait_time <- 1e6
all_stab_results9 <- all_stab_results_small %>%
  distinct() %>%  #there are some duplicated rows....? yes, some treatment combs were repeated when sub1 was made
  filter(num_strains == 9,
         Species == "SB_tot",
         waittime == wait_time) 

all_stab_results9 <- all_stab_results9[with(all_stab_results9, order(var_treat, var_gmax)),]

all_stab_results9$method <- "Simulated"


### Things used for plotting ----



### Fig 5 - log ----

agg_stab_strain9 <- combine_fct("EF_log")

agg_stab_strain9_temp <- agg_stab_strain9 %>%
  dplyr::filter(var_treat=="CB-SB-PB")

agg_stab_strain9 <- agg_stab_strain9 %>%
  dplyr::filter(var_treat!="CB-SB-PB")

agg_stab_strain9 <- rbind(agg_stab_strain9,
                          agg_stab_strain9_temp %>% filter(method %in% c("CB+SBPB","Simulated")) %>% mutate(var_treat="CB-SB-PB & CB+SBPB"),
                          agg_stab_strain9_temp %>% filter(method %in% c("PB+CBSB","Simulated")) %>% mutate(var_treat="CB-SB-PB & PB+CBSB"),
                          agg_stab_strain9_temp %>% filter(method %in% c("SB+CBPB","Simulated")) %>% mutate(var_treat="CB-SB-PB & SB+CBPB"))

agg_stab_strain9 <- agg_stab_strain9 %>%
  mutate(method2 = ifelse(agg_stab_strain9$method=="Simulated","Simulated","Calculated")) %>%
  rename(Calculation=method, Method=method2)%>%
  mutate(Calculation = factor(Calculation, levels = c("Simulated","CB+SB","CB+PB","SB+PB",
                                                      "CB+SBPB","PB+CBSB","SB+CBPB")),
         Method = factor(Method, levels = c("Simulated","Calculated")))



# the legend
fig5legend <- agg_stab_strain9 %>%
  filter(Species == "SB_tot", num_strains==9) %>%
  ggplot(aes(x = stand_var, linetype=Method, col=Calculation)) +
  geom_line(aes(y = hyst_min_log), size = 1) +
  geom_line(aes(y = hyst_max_log), size = 1) +
  scale_color_manual(values = c("black","#e41a1c","#377eb8","#4daf4a","#984ea3","#ff7f00","#a65628")) +
  scale_linetype_manual(values = c(1,2))+
  coord_flip() +
  theme_bw() +
  theme(legend.position="right")

fig5legend <- ggpubr::get_legend(fig5legend)

# the y and x labels
# ylab <- ggplot(data.frame(l = "            Standardised amount of trait variation", x = 1, y = 4.1)) +
ylab <- ggplot(data.frame(l = "Standardised amount of trait variation", x = 1, y = 4.1)) +
  geom_text(aes(x, y, label = l), angle = 90, size=4.5) + 
  theme_void() +
  coord_cartesian(clip = "off")

xlab <- ggplot(data.frame(x = 1, y = 1)) +
  geom_text(aes(x, y), label = expression('log'[10]*"(oxygen diffusivity) (h"^{-1}*")"), size=4.5) + 
  theme_void() +
  coord_cartesian(clip = "off")

# colours and tags used
colours <- c(Simulated="black",`CB+SB`="#e41a1c",`CB+PB`="#377eb8",
             `SB+PB`="#4daf4a",`CB+SBPB`="#984ea3",`PB+CBSB`="#ff7f00",
             `SB+CBPB`="#a65628")

tags <- c(`CB`="a",`SB`="b",`PB`="c",
          `CB-SB`="d",`CB-PB`="e",`SB-PB`="f",
          `CB-SB-PB & CB+SBPB`="g",`CB-SB-PB & PB+CBSB`="h",`CB-SB-PB & SB+CBPB`="i")



# The figure

fig5_log <- split(agg_stab_strain9, f=agg_stab_strain9$var_treat, drop = T)



```

(ref:fig-add-linear) Assessing the additivity or non-additivity of the effects of trait variation. Solid black lines depict the position of the two tipping points at different levels of trait variation, as simulated with the model. The colored, dashed lines show the expected position of tipping points if effects of diversity in different functional groups were additive. The expected tipping point was calculated by adding the effect sizes of trait variation of different functional groups on a linear scale. For example, in (d) the effect size of variation in only CB was added to the effect size of variation in only SB. In (g), the effect size of variation in only CB was added to the effect size of simultaneous variation in both SB and PB.  (CB: cyanobacteria, SB: sulfate-reducing bacteria, PB: phototrophic sulfur bacteria). The dotted vertical lines at -8 and 0 on the x-axis show the extent of the oxygen diffusivity gradient investigated in the simulations.

```{r fig-add-linear, echo = FALSE,  include=TRUE, fig.height=8, fig.width=8, fig.align="center", fig.cap='(ref:fig-add-linear)'}

### Fig 5 - linear ----

agg_stab_strain9 <- combine_fct("EF_lin")

agg_stab_strain9_temp <- agg_stab_strain9 %>%
  dplyr::filter(var_treat=="CB-SB-PB")

agg_stab_strain9 <- agg_stab_strain9 %>%
  dplyr::filter(var_treat!="CB-SB-PB")

agg_stab_strain9 <- rbind(agg_stab_strain9,
                          agg_stab_strain9_temp %>% filter(method %in% c("CB+SBPB","Simulated")) %>% mutate(var_treat="CB-SB-PB & CB+SBPB"),
                          agg_stab_strain9_temp %>% filter(method %in% c("PB+CBSB","Simulated")) %>% mutate(var_treat="CB-SB-PB & PB+CBSB"),
                          agg_stab_strain9_temp %>% filter(method %in% c("SB+CBPB","Simulated")) %>% mutate(var_treat="CB-SB-PB & SB+CBPB"))

agg_stab_strain9 <- agg_stab_strain9 %>%
  mutate(method2 = ifelse(agg_stab_strain9$method=="Simulated","Simulated","Calculated")) %>%
  rename(Calculation=method, Method=method2)%>%
  mutate(Calculation = factor(Calculation, levels = c("Simulated","CB+SB","CB+PB","SB+PB",
                                                      "CB+SBPB","PB+CBSB","SB+CBPB")),
         Method = factor(Method, levels = c("Simulated","Calculated")))



fig5_lin <- split(agg_stab_strain9, f=agg_stab_strain9$var_treat, drop = T)

fig5_lin <- lapply(fig5_lin, function(df){
  tag <- tags[as.character(unique(df$var_treat))]
  plot <- df %>%
    filter(Species == "SB_tot", num_strains==9) %>%
    ggplot(aes(x = stand_var, linetype=Method, col=Calculation)) +
    geom_hline(yintercept = c(-8, 0), col = "black", linetype=3) +
    geom_line(aes(y = hyst_min_log), size = 1) +
    geom_line(aes(y = hyst_max_log), size = 1) +
    facet_wrap( ~ var_treat, scales = "free_y", nrow = 3) +
    labs(y="Oxygen diffusivity\n[log10 uM per hour]",
         x="Amount of trait variation\n[see text for units]",
         linetype="Method",col="Calculation", tag=tag) +
    scale_color_manual(values = colours[unique(df$Calculation)]) +
    scale_linetype_manual(values = c(1,2))+
    coord_flip() +
    theme_bw()+
    theme(legend.position="none",
          axis.title = element_blank(),
          axis.ticks = element_blank(),
          axis.text = element_blank())
  
  if(tag %in% c("a","d")){
    plot <- plot + theme(axis.text.y = element_text(),
                         axis.ticks.y = element_line())
  }
  if(tag %in% c("h","i")){
    plot <- plot + theme(axis.text.x = element_text(),
                         axis.ticks.x = element_line())
  }
  if(tag == "g"){
    plot <- plot + theme(axis.text = element_text(),
                         axis.ticks = element_line())
  }
  plot
})

ylab + Reduce("+",fig5_lin) + fig5legend +
  plot_spacer() + xlab + plot_spacer() +
  plot_layout(ncol = 3, heights = c(25,1), widths = c(1,21,3))


# ggsave(here("experiments/1_figures_main_ms/Figure_5_lin.pdf"),
#        width = 8, height = 8)
```










```{r comp-comb, echo = FALSE, eval = FALSE, include=TRUE, fig.height=8, fig.width=8, fig.align="center", fig.cap='(ref:comp-comb)'}

## This code chunck makes a figure that we used but think does not need
## to be included in the supplement.
## (ref:comp-comb) Another method for comparing the effects of different diversity treatments. Some lines (treatment combinations) are over-plotted and therefore show identical effects of some combinations.


all_stab <- readRDS(here("data",
                         microxanox_version_path,
                         "temporal_method",
                         "event_definition_2",
                         "processed_data",
                         "all_stab_data.RDS"))
wait_time <- 1e6
which_strain <- 9
plot_here <- all_stab %>%
  filter(waittime == wait_time)

  p1 <- plot_here %>%
    #filter(var_treat == "CB") %>%
    filter(num_strains == which_strain) %>%
    filter(Species == "SB_tot") %>%
    ggplot(aes(x = stand_var
               #col = as.factor(num_strains))
    )) +
    geom_line(aes(y = hyst_min_log, col=var_treat), lwd = 1, alpha = 0.8) +
    geom_line(aes(y = hyst_max_log, col=var_treat), lwd = 1, alpha = 0.8) +
    #geom_ribbon(aes(ymin = hyst_min_log, ymax = hyst_max_log),
    #            fill = "green", alpha = 0.1) +
    #facet_grid(var_treat ~ ., scales = "fixed") +
    xlab("Standardised amount of trait variation") +
    ylab(expression(log[10](oxygen~diffusivity)~(h^{-1}))) +
    #labs(fill = "Variation in\nonly these\nfunctional groups") +
    coord_flip() +
    #guides(col = guide_legend(title="Number of strains"),
    #       fill = guide_none()) +
    theme_bw() +
    theme(
      legend.position="top",
      panel.background = element_blank(),
      #strip.text.x = element_blank()
    ) +
    geom_hline(yintercept = c(-8, 0), col = "grey", lwd = 3) #+
    #geom_text(aes(x = 0.95, y = -8.1, label = label), colour = "black") +
    #ggtitle(figure_title)

  p1
 
```

# Log and non-log oxygen diffusivity

For comparability with the [Bush et al publication (2017)](https://www.nature.com/articles/s41467-017-00912-x), we use a log-scale of the oxygen diffusivity gradient when displaying the effect of trait variation on tipping point positions. In our work, we only make comparisons where the outcome is the same when a log- or an arithmetic-scale of oxygen diffusivity was used. For example, we compare the direction of the diversity effects of the three functional groups on the resilience of the ecosystem state that they dominate (CB: oxic state, i.e. oxic-anoxic tipping point; SB and PB: anoxic state, i.e. anoxic-oxic tipping point). The direction of the effects does not depend on the log/arithmetic scale decision.

In contrast, we avoid comparisons where the outcome would differ when a log- or an arithmetic-scale of oxygen diffusivity was used. An example of a comparison that we avoid is of the magnitude of the diversity effect between the oxic-anoxic transition and the anoxic-oxic transition. We avoid this comparison of magnitudes because with a log scale the diversity effect is larger for the oxic-anoxic transition (Figure \@ref(fig:fig-log-nonlog)n), while on an arithmetic scale the diversity effect is larger for the anoxic-oxic transition (Figure \@ref(fig:fig-log-nonlog)o). 


(ref:fig-log-nonlog) The effect of trait diversity (y-axis) viewed with a log axis (left panels) or arithmetic axis (right panels) for the oxygen diffusivity gradient (x-axis). The positions of the tipping points (anoxic to oxic in blue, oxic to anoxic in red) and the extent of the region of bistability (green or orange shaded region) are shown. The left-most region is anoxic only, the middle is oxic or anoxic depending on historical conditions, and the right-most is oxic only. In the left panels the thick grey vertical lines at -8 and 0 on the x-axis show the extent of the oxygen diffusivity gradient investigated in the simulations. In panels e, f, i, and j the orange region shows when patterns of coexistence were atypical of the classical pattern shown in Figure 1a of the main text (e.g. when cyanobacteria coexisted with sulfur bacteria at low oxygen diffusivities). In the right-hand panels (arithmetic x-axis scale) the blue anoxic to oxic tipping point line often has steps. This is because the simulated oxygen diffusivity levels were evenly distributed on the log-oxygen diffusivity scale, and so are relatively far apart at high oxygen diffusivities on the arithmetic scale. 

```{r fig-log-nonlog, echo = FALSE,  include=TRUE, fig.height=8, fig.width=8, fig.align="center", fig.cap='(ref:fig-log-nonlog)'}
all_stab <- readRDS(here("data",
                         microxanox_version_path,
                         "temporal_method",
                         "event_definition_2",
                         "processed_data",
                         "all_stab_data.RDS"))

wait_time <- 1e6
num_strains <- 9
plot_here <- all_stab %>%
  filter(waittime == wait_time)

plot_here <- plot_here %>%
  mutate(stab_states_type = "standard",
         stab_states_type = case_when(
           (var_treat == "CB" & stand_var > 1.00) ~ "non-standard",
           (var_treat == "SB" & stand_var > 1.00) ~ "non-standard",
           (var_treat == "PB" & stand_var > 0.65) ~ "non-standard",
           (var_treat == "CB-SB" & stand_var > 1.00) ~ "non-standard",
           (var_treat == "CB-PB" & stand_var > 0.65) ~ "non-standard",
           (var_treat == "SB-PB" & stand_var > 1.00) ~ "non-standard",
           (var_treat == "CB-SB-PB" & stand_var > 1.00) ~ "non-standard"
         ),
         stab_states_type = ifelse(is.na(stab_states_type), "standard", "non-standard")
  )


p1 <- fig_div_vs_o2diff_1strain_7row_stable_state_type(plot_here,
                                                       which_strain = num_strains,
                                                       figure_title = NULL
                                                       # figure_title = paste0(wait_time,
                                                       #                       " wait time,\n",
                                                       #                       num_strains,
                                                       #                       " strains")
)


p2 <- fig_div_vs_o2diff_1strain_7row_stable_state_type_nonlog(plot_here,
                                                       which_strain = num_strains,
                                                       figure_title = NULL
                                                       # figure_title = paste0(wait_time,
                                                       #                       " wait time,\n",
                                                       #                       num_strains,
                                                       #                       " strains")
)

p1 + p2


```


# References

Bush, T., Diao, M., Allen, R.J., Sinnige, R., Muyzer, G. & Huisman, J. (2017). Oxic-anoxic regime shifts mediated by feedbacks between biogeochemical processes and microbial community dynamics. Nat. Commun., 8, 789.

Cardinale, B.J., Duffy, J.E., Gonzalez, A., Hooper, D.U., Perrings, C., Venail, P., et al. (2012). Biodiversity loss and its impact on humanity. Nature, 486, 59–67.

Cardinale, B.J., Matulich, K.L., Hooper, D.U., Byrnes, J.E., Duffy, E., Gamfeldt, L., et al. (2011). The functional role of producer diversity in ecosystems. Am. J. Bot., 98, 572–592.

Cardinale, B.J., Srivastava, D.S., Emmett Duffy, J., Wright, J.P., Downing, A.L., Sankaran, M., et al. (2006). Effects of biodiversity on the functioning of trophic groups and ecosystems. Nature, 443, 989–992.

Fargione, J., Tilman, D., Dybzinski, R., Lambers, J.H.R., Clark, C., Harpole, W.S., et al. (2007). From selection to complementarity: shifts in the causes of biodiversity–productivity relationships in a long-term biodiversity experiment. Proc. R. Soc. B Biol. Sci., 274, 871–876.

Hughes, A.R. & Stachowicz, J.J. (2011). Seagrass genotypic diversity increases disturbance response via complementarity and dominance. J. Ecol., 99, 445–453.

Loreau, M., Barbier, M., Filotas, E., Gravel, D., Isbell, F., Miller, S.J., et al. (2021). Biodiversity as insurance: from concept to measurement and application. Biol. Rev., 96, 2333–2354.

Loreau, M. & Hector, A. (2001). Partitioning selection and complementarity in biodiversity experiments. Nature, 412, 72–76.

Loreau, M. & Mazancourt, C. de. (2013). Biodiversity and ecosystem stability: a synthesis of underlying mechanisms. Ecol. Lett., 16, 106–115.

Luther, G.W., Findlay, A., MacDonald, D., Owings, S., Hanson, T., Beinart, R., et al. (2011). Thermodynamics and kinetics of sulfide oxidation by oxygen: a look at inorganically controlled reactions and biologically mediated processes in the environment. Front. Microbiol., 2, 62.

Moisset de Espanés, P., Ramos-Jiliberto, R. & Soto, J.A. (2021). Avoiding artifacts when varying the number of species in ecological models. Ecol. Lett., 24, 1976–1987.

Steiner, C.F., Long, Z.T., Krumins, J.A. & Morin, P.J. (2006). Population and Community Resilience in Multitrophic Communities. Ecology, 87, 996–1007.


