User guide – microxanox R Package for microbial oxic and anoxic ecosystem simulations
Package version
The microxanox package may be update / changed, so please take care to work with a specific version. This user guide is requires at least version 0.4, and the loaded version is:
## [1] '0.9.3'Introduction
The aims of the microxanox package are:
- Facilitate reproduction of the results of Bush et al 2017.
- Include new functionality that allows new research questions to be answered.
We added functionality for:
- Adding random noise in substrate concentrations.
- Adding temporally varying oxygen diffusivity.
- Including immigration.
- Setting minimum population abundances.
- Multiple strains per functional group.
We also include a function for visualising results. It is provided for convenience, and not for production of perfectly formatted print-ready graphs.
Here (again) is the original publication Bush et al (2017): https://www.nature.com/articles/s41467-017-00912-x
And here is the Supplementary Information, including the table of
parameter values as implemented here when using the bush
parameter: https://static-content.springer.com/esm/art%3A10.1038%2Fs41467-017-00912-x/MediaObjects/41467_2017_912_MOESM1_ESM.pdf
Notation
Note that in this package we use the notation from the Table S1 of
the Bush
et al 2017 Supplementary information, and not the notation in the
main text of Bush et al
2017. (There is only one small difference between these: In the main
text the notation is simplified to y_S_PB (yS, PB)
(which is y_SR_PB (ySR, PB)
in Table S1) and y_S_SB (yS, SB)
(which is y_SO_SB (ySO, SB)
in Table S1). We also use Pr_cB (PrB)
instead of P_CB (PCB).
Model dimensions
- Time: hours
- Volume: litres
- Substrate quantity: micromoles
- Organism quantity: cells
Modelling framework
The framework used for running a simulation, and for varying aspects of a simulation run is:
- Define a parameter set of class
runsim_parameter(i.e. parameters and initial conditions) by using the functionsnew_runsim_parameter(),new_initial_state()andnew_strain_parameter(). - The function call
run_simulation(parameter)would therefore run a simulation using the parameter as defined in the objectparameter.
Conditions that must be defined
All parameter needed to run the simulations are encapsulated in an S3
object of class runsim_parameter. This object is returned
by the function new_runsim_parameter().
## $dynamic_model
## [1] NA
##
## $event_definition
## [1] NA
##
## $strain_parameter
## [1] NA
##
## $event_interval
## [1] NA
##
## $noise_sigma
## [1] NA
##
## $minimum_abundances
## [1] NA
##
## $sim_duration
## [1] NA
##
## $sim_sample_interval
## [1] NA
##
## $log10a_series
## [1] NA
##
## $asym_factor
## [1] FALSE
##
## $solver_method
## [1] "radau"
##
## attr(,"class")
## [1] "list" "runsim_parameter"All parameter values are set initially to NA (except the solver
method, which is set to radau) and therefore need to be set. This can
either be done by providing the parameters to be set as named arguments
to the function new_runsim_parameter():
## $dynamic_model
## function (t, state, parameters, log10a_forcing_func, ...)
## {
## CB <- state[grep("CB", names(state))]
## names_CB <- names(CB)[order(names(CB))]
## CB <- as.numeric(CB[order(names(CB))])
## PB <- state[grep("PB", names(state))]
## names_PB <- names(PB)[order(names(PB))]
## PB <- as.numeric(PB[order(names(PB))])
## SB <- state[grep("SB", names(state))]
## names_SB <- names(SB)[order(names(SB))]
## SB <- as.numeric(SB[order(names(SB))])
## CB_growth_rate <- growth1(state["P"], parameters$CB$g_max_CB,
## parameters$CB$k_CB_P) * inhibition(state["SR"], parameters$CB$h_SR_CB) *
## CB
## CB_mortality_rate <- parameters$CB$m_CB * CB
## CB_rate <- CB_growth_rate - CB_mortality_rate + parameters$CB$i_CB
## PB_growth_rate <- growth2(state["P"], state["SR"], parameters$PB$g_max_PB,
## parameters$PB$k_PB_P, parameters$PB$k_PB_SR) * inhibition(state["O"],
## parameters$PB$h_O_PB) * PB
## PB_mortality_rate <- parameters$PB$m_PB * PB
## PB_rate <- PB_growth_rate - PB_mortality_rate + parameters$PB$i_PB
## SB_growth_rate <- growth2(state["P"], state["SO"], parameters$SB$g_max_SB,
## parameters$SB$k_SB_P, parameters$SB$k_SB_SO) * inhibition(state["O"],
## parameters$SB$h_O_SB) * SB
## SB_mortality_rate <- parameters$SB$m_SB * SB
## SB_rate <- SB_growth_rate - SB_mortality_rate + parameters$SB$i_SB
## SO_rate <- sum(1/parameters$PB$y_SR_PB * PB_growth_rate) -
## sum(1/parameters$SB$y_SO_SB * SB_growth_rate) + parameters$c *
## state["O"] * state["SR"] + parameters$a_S * (parameters$back_SO -
## state[["SO"]])
## SR_rate <- -sum(1/parameters$PB$y_SR_PB * PB_growth_rate) +
## sum(1/parameters$SB$y_SO_SB * SB_growth_rate) - parameters$c *
## state["O"] * state["SR"] + parameters$a_S * (parameters$back_SR -
## state["SR"])
## O_rate <- sum(parameters$CB$Pr_CB * CB_growth_rate) - parameters$c *
## state["O"] * state["SR"] + 10^log10a_forcing_func(t) *
## (parameters$back_O - state["O"])
## P_rate <- -sum(1/parameters$CB$y_P_CB * CB_growth_rate) -
## sum(1/parameters$PB$y_P_PB * PB_growth_rate) - sum(1/parameters$SB$y_P_SB *
## SB_growth_rate) + parameters$a_P * (parameters$back_P -
## state["P"])
## result <- list(c(CB_rate, PB_rate, SB_rate, SO_rate = SO_rate,
## SR_rate = SR_rate, O_rate = O_rate, P_rate = P_rate),
## a = log10a_forcing_func(t))
## names(result[[1]]) <- c(parameters$CB$strain_name, parameters$PB$strain_name,
## parameters$SB$strain_name, "SO_rate", "SR_rate", "O_rate",
## "P_rate")
## return(result)
## }
## <bytecode: 0x56413d3a7138>
## <environment: namespace:microxanox>
##
## $event_definition
## [1] NA
##
## $strain_parameter
## [1] NA
##
## $event_interval
## [1] NA
##
## $noise_sigma
## [1] NA
##
## $minimum_abundances
## [1] NA
##
## $sim_duration
## [1] NA
##
## $sim_sample_interval
## [1] NA
##
## $log10a_series
## [1] NA
##
## $asym_factor
## [1] FALSE
##
## $solver_method
## [1] "radau"
##
## attr(,"class")
## [1] "list" "runsim_parameter"Parameters can also be set after creation of a parameter set:
This parameter set has the following elements:
## $dynamic_model
## function (t, state, parameters, log10a_forcing_func, ...)
## {
## CB <- state[grep("CB", names(state))]
## names_CB <- names(CB)[order(names(CB))]
## CB <- as.numeric(CB[order(names(CB))])
## PB <- state[grep("PB", names(state))]
## names_PB <- names(PB)[order(names(PB))]
## PB <- as.numeric(PB[order(names(PB))])
## SB <- state[grep("SB", names(state))]
## names_SB <- names(SB)[order(names(SB))]
## SB <- as.numeric(SB[order(names(SB))])
## CB_growth_rate <- growth1(state["P"], parameters$CB$g_max_CB,
## parameters$CB$k_CB_P) * inhibition(state["SR"], parameters$CB$h_SR_CB) *
## CB
## CB_mortality_rate <- parameters$CB$m_CB * CB
## CB_rate <- CB_growth_rate - CB_mortality_rate + parameters$CB$i_CB
## PB_growth_rate <- growth2(state["P"], state["SR"], parameters$PB$g_max_PB,
## parameters$PB$k_PB_P, parameters$PB$k_PB_SR) * inhibition(state["O"],
## parameters$PB$h_O_PB) * PB
## PB_mortality_rate <- parameters$PB$m_PB * PB
## PB_rate <- PB_growth_rate - PB_mortality_rate + parameters$PB$i_PB
## SB_growth_rate <- growth2(state["P"], state["SO"], parameters$SB$g_max_SB,
## parameters$SB$k_SB_P, parameters$SB$k_SB_SO) * inhibition(state["O"],
## parameters$SB$h_O_SB) * SB
## SB_mortality_rate <- parameters$SB$m_SB * SB
## SB_rate <- SB_growth_rate - SB_mortality_rate + parameters$SB$i_SB
## SO_rate <- sum(1/parameters$PB$y_SR_PB * PB_growth_rate) -
## sum(1/parameters$SB$y_SO_SB * SB_growth_rate) + parameters$c *
## state["O"] * state["SR"] + parameters$a_S * (parameters$back_SO -
## state[["SO"]])
## SR_rate <- -sum(1/parameters$PB$y_SR_PB * PB_growth_rate) +
## sum(1/parameters$SB$y_SO_SB * SB_growth_rate) - parameters$c *
## state["O"] * state["SR"] + parameters$a_S * (parameters$back_SR -
## state["SR"])
## O_rate <- sum(parameters$CB$Pr_CB * CB_growth_rate) - parameters$c *
## state["O"] * state["SR"] + 10^log10a_forcing_func(t) *
## (parameters$back_O - state["O"])
## P_rate <- -sum(1/parameters$CB$y_P_CB * CB_growth_rate) -
## sum(1/parameters$PB$y_P_PB * PB_growth_rate) - sum(1/parameters$SB$y_P_SB *
## SB_growth_rate) + parameters$a_P * (parameters$back_P -
## state["P"])
## result <- list(c(CB_rate, PB_rate, SB_rate, SO_rate = SO_rate,
## SR_rate = SR_rate, O_rate = O_rate, P_rate = P_rate),
## a = log10a_forcing_func(t))
## names(result[[1]]) <- c(parameters$CB$strain_name, parameters$PB$strain_name,
## parameters$SB$strain_name, "SO_rate", "SR_rate", "O_rate",
## "P_rate")
## return(result)
## }
## <bytecode: 0x56413d3a7138>
## <environment: namespace:microxanox>
##
## $event_definition
## [1] NA
##
## $strain_parameter
## [1] NA
##
## $event_interval
## [1] NA
##
## $noise_sigma
## [1] NA
##
## $minimum_abundances
## [1] NA
##
## $sim_duration
## [1] NA
##
## $sim_sample_interval
## [1] NA
##
## $log10a_series
## [1] NA
##
## $asym_factor
## [1] FALSE
##
## $solver_method
## [1] "radau"
##
## attr(,"class")
## [1] "list" "runsim_parameter"Now we will discuss elements in the parameter set that must be set, with an example of setting them:
The dynamic model
The ordinary differential equations for the rates of change are
specified in the function bushplus_dynamic_model(). This
augmented version of the model published in [Bush et al
2017](https://www.nature.com/articles/s41467-017-00912-x) can handle
multiple strains within each of the three functional groups, temporal
variation in oxygen diffusivity, and events (see below for event
definition).
The parameter and initial state values
The parameter values and initial state values are held in one object.
The function new_starter can be used to create this object.
In the follow, we specify one strain per functional group, all parameter
values as defined in [Bush et al
2017](https://www.nature.com/articles/s41467-017-00912-x), and
starting conditions as the anoxic favouring conditions in Bush et al
2017.
parameter_set1$strain_parameter <- new_strain_parameter(
n_CB = 1,
values_CB = "bush",
n_PB = 1,
values_PB = "bush",
n_SB = 1,
values_SB = "bush",
values_other = "bush",
values_initial_state = "bush_anoxic_fig2ab"
)Here, n_CB, n_{PB}, and n_{SB}
are the number of strains within each functional group and
values_{CB} (etc) are the parameter values of the strains,
i.e. "bush" uses the parameter as used in the Bush paper.
The same applies to values_initial_state, which sets the
initial state of the system to a preset from the Bush et al
2017 publication when using "bush_anoxic_fig2ab".
Here are the parameter values of the CB strain, for example
## strain_name g_max_CB k_CB_P h_SR_CB y_P_CB Pr_CB m_CB i_CB
## 1 CB_1 0.05 0.2 300 1.67e+08 6e-09 0.02 0and these are initial state values used to start a simulation:
## CB_1 PB_1 SB_1 SO SR O P
## 5e+01 1e+07 1e+07 3e+02 3e+02 1e+01 1e+01
## attr(,"class")
## [1] "initial_state" "numeric"Modifying initial states and parameter values
Note that the function new_initial_state can be used to
get a couple of preset initial states:
## CB_1 PB_1 SB_1 SO SR O P
## 1e+08 1e+02 1e+02 5e+02 5e+01 3e+02 4e+00
## attr(,"class")
## [1] "initial_state" "numeric"And the function new_CB_strain_parameter (and SB and PB
equivalents) can be used to get preset organism parameters. Currently
there is only one preset, however: “bush”.
## strain_name g_max_CB k_CB_P h_SR_CB y_P_CB Pr_CB m_CB i_CB
## 1 CB_1 0.05 0.2 300 1.67e+08 6e-09 0.02 0The sampling interval
How frequently are the values of state variables recorded:
The oxygen diffusivity dynamics
Specified by giving a vector of values of aO. These are then automatically spaced equally along the duration of the simulation. Here we specify only two values of aO, and so these are the start and end values, and they are the same. Hence there is no temporal change in aO.
parameter_set1$log10a_series <- c(
log10(parameter_set1$strain_parameter$a_O),
log10(parameter_set1$strain_parameter$a_O)
)Note that here we use the already set default value of aO and therefore must have already defined the default parameter values.
Events
These are “things” that happen during the simulation. One has to specify what happens, and when it happens.
What happens
There are various available event definitions. Please see the functions for more information. Even definition 1 can add noise to a state variable, and/or assure that abundances are limited to a minimum:
Default noise event options
Default to no noise added to substrate concentrations:
Setting parameters
The parameters can also be set in the command that creates the
run_sim_parameter object:
sp <- new_strain_parameter(
n_CB = 1,
values_CB = "bush",
n_PB = 1,
values_PB = "bush",
n_SB = 1,
values_SB = "bush",
values_other = "bush",
values_initial_state = "bush_anoxic_fig2ab"
)
parameter_set1 <- new_runsim_parameter(
dynamic_model = bushplus_dynamic_model,
strain_parameter = sp,
sim_duration = 3600,
sim_sample_interval = 1,
log10a_series = c(
log10(sp$a_O),
log10(sp$a_O)
),
event_definition = event_definition_1,
event_interval = 10,
noise_sigma = 0,
minimum_abundances = rep(0, 3)
)
names(parameter_set1$minimum_abundances) <- c("CB", "PB", "SB")
rm(sp)The advantage of this approach is that the names of the parameter are checked and an error is raised when it does not match, e.g.
## Error in new_runsim_parameter(dynamical_model = bushplus_dynamic_model): Non defined parameter supplied: dynamical_modelRun some example simulations
1 strain anoxic start
Once all the conditions are assigned to the parameter object, then
these can be given to the run_simulation():
This will return an object of class run_sim_results
which is identical to parameter, except that it contains
one additional element, namely result, which contains the
results of the simulation.
Consequently, saving anoxic_start_sim_results will save
all parameters used to (re)run the simulation as well as the result:
As a sidenote, the simulation can be run again by simply calling
The convenience function plot_dynamics can be used to
easily visualise the dynamics of a simulation:
1 strain oxic start
oxic_start_sim_results <- run_simulation(parameter_set1)
saveRDS(oxic_start_sim_results, file = mdat("oxic_start_sim_results.RDS"))
Three strains per functional group
Here we simulate a system with three strains per functional group. The system dynamics should not (and do not) differ from the 1 strain version, because the strains within a functional group are (in this case) identical.
For example, here are the parameters of the three CB strains:
## strain_name g_max_CB k_CB_P h_SR_CB y_P_CB Pr_CB m_CB i_CB
## 1 CB_1 0.05 0.2 300 1.67e+08 6e-09 0.02 0
## 2 CB_2 0.05 0.2 300 1.67e+08 6e-09 0.02 0
## 3 CB_3 0.05 0.2 300 1.67e+08 6e-09 0.02 0And the initial state values
## CB_1 CB_2 CB_3 PB_1 PB_2 PB_3 SB_1 SB_2 SB_3 SO SR O P
## 5e+01 5e+01 5e+01 1e+07 1e+07 1e+07 1e+07 1e+07 1e+07 3e+02 3e+02 1e+01 1e+01
## attr(,"class")
## [1] "initial_state" "numeric"So we can see the different strains, we now give them slightly different initial abundances:
parameter_set1$strain_parameter$initial_state["CB_1"] <- 5e3
parameter_set1$strain_parameter$initial_state["CB_2"] <- 1e3
parameter_set1$strain_parameter$initial_state["CB_3"] <- 0
parameter_set1$strain_parameter$initial_state["PB_1"] <- 1e7
parameter_set1$strain_parameter$initial_state["PB_2"] <- 1e6
parameter_set1$strain_parameter$initial_state["PB_3"] <- 1e5
parameter_set1$strain_parameter$initial_state["SB_1"] <- 1e7
parameter_set1$strain_parameter$initial_state["SB_2"] <- 1e6
parameter_set1$strain_parameter$initial_state["SB_3"] <- 1e5
# parameter$strain_parameter$initial_state["SO"] <- 300
# parameter$strain_parameter$initial_state["SR"] <- 300
# parameter$strain_parameter$initial_state["O"] <- 1e1
# parameter$strain_parameter$initial_state["P"] <- 1e1parameter_set1$initial_state <- parameter_set1$strain_parameter$initial_state
simulation_result3 <- run_simulation(parameter_set1)
saveRDS(simulation_result3, file = mdat("multistrain3.RDS"))
Get steady state(s)
There are two methods for finding steady states. The first runs a separate simulation for each combination of starting conditions and oxygen diffusivity (let us term this the Replication method). The second runs only two simulations, with step-wise and slowly temporally increasing or decreasing oxygen diffusivities and recorded of state just before change to a new oxygen diffusivity (let us term this the Temporal method).
Replication method.
1 strain per functional group
Finding the steady state of the system is a common task. The function
run_replication_ssfind() finds the state of a system after
a certain length of simulation, for different values of the a0 (oxygen diffusivity)
parameter, and for different initial conditions (N).
Note: The function run_replication_ssfind() can be told
to use multiple cores by setting the relevant option, e.g.,
options(mc.cores=8)
We first set the conditions that will apply to all simulations. Note
that these are automatically picked up by the functions called by
run_replication_ssfind() and so do not need to be
passed.
We will work with a one strain per functional group model, and unless otherwise changed (which we do) starting with oxic favouring initial conditions.
ma <- rep(0, 3)
names(ma) <- c("CB", "PB", "SB")
parameter_set2 <- new_replication_ssfind_parameter(
dynamic_model = bushplus_dynamic_model,
strain_parameter = new_strain_parameter(
n_CB = 1,
n_PB = 1,
n_SB = 1,
values_initial_state = "bush_ssfig3"
),
sim_duration = 50000,
sim_sample_interval = 50000,
event_definition = event_definition_1,
event_interval = 50000,
ss_expt = NA,
noise_sigma = 0,
minimum_abundances = ma
)
rm(ma)Next we need to define the values of a0 and initial conditions that we will get the steady states of. We create quite continuous variation a0. We, however, only need to vary the initial abundance of CB and have it either high or low, in order to find which of the two possible states that can occur do occur.
grid_num_a <- 100 ## number of a_0 values
grid_num_N <- 2 ## number of N values
initial_CBs <- 10^seq(0, 10, length=grid_num_N) ## sequence of N values
initial_PBs <- 1e8 ## not varied
initial_SBs <- 1e8 ## not varied
a_Os <- 10^seq(-6, -2, length=grid_num_a) ## sequence of a_0 values
## next line creates all possible combinations
expt <- expand.grid(N_CB = initial_CBs,
N_PB = initial_PBs,
N_SB = initial_SBs,
a_O = a_Os)
parameter_set2$ss_expt <- expt
rm(grid_num_a, grid_num_N, initial_CBs, initial_PBs, initial_SBs, a_Os, expt)
# next line runs the "experiment", calling various functions to do this.
sim_res <- run_replication_ssfind(parameter_set2)
saveRDS(sim_res, file = mdat("ss_1strain.RDS"))And the results:
p1 <- ss_1strain$result %>%
mutate(direction = ifelse(initial_N_CB == 1, "up", "down")) %>%
select(a, direction,
starts_with("CB"),
starts_with("PB"),
starts_with("SB")) %>%
gather(key = Species, value=Quantity, 3:ncol(.)) %>%
ggplot() +
geom_path(aes(x = a, y = log10(Quantity+1), col = Species, linetype = direction),
size=1) +
ylab(expression(Log[10](Abundance+1))) +
xlab(expression(Log[10](Oxygen~diffusivity))) +
ylim(0,10) +
#theme(plot.title=element_text(size=10, hjust=0.5)) +
labs(title="Organisms")## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.##ggsave("figures/CB_SB_PB.png", width = 3, height = 2)
p2 <- ss_1strain$result %>%
mutate(direction = ifelse(initial_N_CB == 1, "up", "down")) %>%
select(a, direction ,SO, SR, O, P) %>%
gather(key = Species, value=Quantity, 3:ncol(.)) %>%
ggplot() +
geom_path(aes(x = a, y = log10(Quantity), col = Species, linetype = direction)) +
labs(title="Substrates")
p1 / p2
3 strains per functional group
And now finding the steady states of a three strain per functional group system, though with identical strains in each group. Hence the steady states should be the same as the one strain version.
Note that the run_replication_ssfind() function divides
the specified starting abundance of a functional group and distributes
that equally among the strains of that functional group.
parameter_set2$strain_parameter <- new_strain_parameter(
n_CB = 3,
n_PB = 3,
n_SB = 3,
values_initial_state = "bush_ssfig3"
)
ss_3strain <- run_replication_ssfind(parameter_set2)
saveRDS(ss_3strain, file = mdat("ss_3strain.RDS"))And the results (plotting code hidden this time):
Temporal method
The temporal method involves two simulations for a particular system configuration (parameter set). In one simulation the oxygen diffusivity is increased in a step-wise fashion. In the other it is decreased in a step-wise fashion. That is, oxygen diffusivity is held at a constant level for long enough for steady state to be reach, that state is recorded, and then a slightly higher (or lower) oxygen diffusivity value is set. Hence, at that time point, the system is effectively started with initial conditions that are the state of the system in the previous time step. We hereafter term this “the temporal method”. The results below are for a very short simulation, hence the region of bistability is large. I.e. the rate of temporal change is fast, so the system lags, and this creates a larger region of bistability.
log10a_series <- seq(-8, 0, length = 100)
wait_time <- 1e3 # time spent at each step
num_strains <- 3
total_initial_abundances <- 10^5
num_CB_strains <- num_strains
num_SB_strains <- num_strains
num_PB_strains <- num_strains
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 = 1000,
noise_sigma = 0,
minimum_abundances = rep(1, 3),
strain_parameter = sp,
log10a_series = log10a_series
)
names(parameter$minimum_abundances) <- c("CB", "PB", "SB")
rm(sp)
parameter$sim_duration <- wait_time * length(parameter$log10a_series)
parameter$sim_sample_interval <- wait_time
parameter <- set_temporal_ssfind_initial_state(
parameter,
total_initial_abundances,
total_initial_abundances,
total_initial_abundances
)
## be careful, this could take a long time to run if parameters above are changed (the demonstration parameters are with very short wait time).
temporal_ss_find_result <- run_temporal_ssfind(parameter)
saveRDS(temporal_ss_find_result, file = mdat("temporal_ss_find.RDS"))temporal_ss_find_result <- readRDS(file = mdat("temporal_ss_find.RDS"))
p1 <- temporal_ss_find_result %>%
mutate(a = 10^a_O) %>%
gather(species, quantity, 2:(ncol(.) - 2)) %>%
mutate(log10_quantity = log10(quantity + 1)) %>%
dplyr::filter(species == "O") %>%
ggplot(aes(x = log10(a), y = log10_quantity, col = species)) +
geom_path() +
ylab("Log(Quantity)") +
xlab("Oxygen diffusivity") +
theme(legend.position="none") +
geom_vline(xintercept = c(-8, 0), col = "grey", lwd = 3)## Warning: There was 1 warning in `mutate()`.
## ℹ In argument: `log10_quantity = log10(quantity + 1)`.
## Caused by warning:
## ! NaNs produced
Get measures of nonlinearity and hysteresis
Replication method
The get_stability_measures() function takes an object of
the type returned by the run_replication_ssfind() function
and returns various measures of the stability of each of the state
variables.
The measures of non-linearity and hysteresis were used in Garnier et al (2020) (link to online article).
(Warning: Care is required when processing the produced data. Some log transformations with addition of 1 are used; NAs are produced in some cases; and some measures are affected by the number of strains, even when the strains are identical.)
For example:
## Warning: Returning data frames from `filter()` expressions was deprecated in dplyr
## 1.0.8.
## ℹ Please use `if_any()` or `if_all()` instead.
## ℹ The deprecated feature was likely used in the microxanox package.
## Please report the issue at <https://github.com/UZH-PEG/microxanox/issues>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.## # A tibble: 7 × 13
## Species hyst_tot_raw hyst_range_raw hyst_min_raw hyst_max_raw nl_up_raw
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 CB_1 4.07e+8 0.00154 0.0000123 0.00156 1.55e+7
## 2 O 2.11e+1 0 0 0 4.17e+0
## 3 P 6.06e-2 0 0 0 5.95e-3
## 4 PB_1 4.57e+7 0.00141 0.0000102 0.00142 3.07e+6
## 5 SB_1 1.70e+8 0.00155 0.0000102 0.00156 6.69e+6
## 6 SO 1.34e+2 0 0 0 7.38e+0
## 7 SR 1.34e+2 0 0 0 7.38e+0
## # ℹ 7 more variables: nl_down_raw <dbl>, hyst_tot_log <dbl>,
## # hyst_range_log <dbl>, hyst_min_log <dbl>, hyst_max_log <dbl>,
## # nl_up_log <dbl>, nl_down_log <dbl>And with three strains per functional group:
## # A tibble: 13 × 13
## Species hyst_tot_raw hyst_range_raw hyst_min_raw hyst_max_raw nl_up_raw
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 CB_1 1.36e+8 0.00154 0.0000123 0.00156 5163765.
## 2 CB_2 1.36e+8 0.00154 0.0000123 0.00156 5163765.
## 3 CB_3 1.36e+8 0.00154 0.0000123 0.00156 5163765.
## 4 O 2.11e+1 0 0 0 4.17
## 5 P 6.06e-2 0 0 0 0.00595
## 6 PB_1 1.52e+7 0.00128 0.0000112 0.00129 1023085.
## 7 PB_2 1.52e+7 0.00128 0.0000112 0.00129 1023085.
## 8 PB_3 1.52e+7 0.00128 0.0000112 0.00129 1023085.
## 9 SB_1 5.65e+7 0.00154 0.0000112 0.00156 2229592.
## 10 SB_2 5.65e+7 0.00154 0.0000112 0.00156 2229592.
## 11 SB_3 5.65e+7 0.00154 0.0000112 0.00156 2229592.
## 12 SO 1.34e+2 0 0 0 7.38
## 13 SR 1.34e+2 0 0 0 7.38
## # ℹ 7 more variables: nl_down_raw <dbl>, hyst_tot_log <dbl>,
## # hyst_range_log <dbl>, hyst_min_log <dbl>, hyst_max_log <dbl>,
## # nl_up_log <dbl>, nl_down_log <dbl>Temporal method
All as in the previous section, except that
get_stability_measures() function takes an object of the
type returned by the run_temporal_ssfind function.
## # A tibble: 13 × 13
## Species hyst_tot_raw hyst_range_raw hyst_min_raw hyst_max_raw nl_up_raw
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 CB_1 137364902. 8 -8 0 6503411.
## 2 CB_2 137364902. 8 -8 0 6503411.
## 3 CB_3 137364902. 8 -8 0 6503411.
## 4 O 20.0 0 0 0 6.95
## 5 P 0.102 0 0 0 0.0100
## 6 PB_1 19270729. 5.49 -8 -2.51 802364.
## 7 PB_2 19270729. 5.49 -8 -2.51 802364.
## 8 PB_3 19270729. 5.49 -8 -2.51 802364.
## 9 SB_1 56879841. 5.49 -8 -2.51 2765768.
## 10 SB_2 56879841. 5.49 -8 -2.51 2765768.
## 11 SB_3 56879841. 5.49 -8 -2.51 2765768.
## 12 SO 126. 0 0 0 8.48
## 13 SR 126. 0 0 0 8.48
## # ℹ 7 more variables: nl_down_raw <dbl>, hyst_tot_log <dbl>,
## # hyst_range_log <dbl>, hyst_min_log <dbl>, hyst_max_log <dbl>,
## # nl_up_log <dbl>, nl_down_log <dbl>