,
Rainer M. Krug [ctb] ,
Pascal Baertschi [ctb] | Title: | Oxic-Anoxic Regime Shifts in Microbial Communities |
|---|---|
| Description: | Model to simulate a three functional group system with four chemical substrates using a set of ordinary differential equations. Simulations can be run individually or over a parameter range, to find stable states. The model features multiple species per functional group, where the number is only limited by computational constraints. The R package is constructed in such a way, that the results contain the input parameter used, so that a saved results can be loaded again and thesimulation be repeated. |
| Authors: | Owen L. Petchey [aut, cre] ,
Rainer M. Krug [ctb] ,
Pascal Baertschi [ctb] |
| Maintainer: | Owen L. Petchey <[email protected]> |
| License: | MIT + file LICENSE |
| Version: | 0.9.3 |
| Built: | 2024-11-11 06:02:21 UTC |
| Source: | https://github.com/UZH-PEG/microxanox |
Ecosystems containing microbes can be oxic (oxygen present) or anoxic (oxygen absent). Under some conditions these states can be both stable even with identical biotic and abiotic circumstances, i.e. they can be alternate stable states, implying history alone can determine the present state. The biological reason is that the microbes can alter the environmental conditions in a way that favour themselves, while making them less suitable for other organisms. This creates mutual inhibition. Which state occurs also depends on numerous abiotic and biotic factors. This package facilitates simulation studies of the influence of several of these factors, and is based on work described in (Bush et al) 2017 (DOI: 10.1038/s41467-017-00912-x)
Vignettes include:
- A user guide. - A reproduction of some results of Bush et al 2017, on which the simulation is based. It shows the alternate stable states. It also extends the original study by adding temporal forcing and therefore temporal switching between states, and also effects of biotic composition on the response to environmental change.
Add variability as defined in the argument variability to the variables as follows:
strain_parameter$CB$g_max_CB <- variability(strain_parameter$CB$g_max_CB, CB_var_gmax)
strain_parameter$CB$h_SR_CB <- variability(strain_parameter$CB$h_SR_CB, CB_var_h )
strain_parameter$SB$g_max_SB <- variability(strain_parameter$SB$g_max_SB, SB_var_gmax)
strain_parameter$SB$h_O_SB <- variability(strain_parameter$SB$h_O_SB, SB_var_h )
strain_parameter$PB$g_max_PB <- variability(strain_parameter$PB$g_max_PB, PB_var_gmax)
strain_parameter$PB$h_O_PB <- variability(strain_parameter$PB$h_O_PB, PB_var_h )
add_strain_var( strain_parameter, CB_var_gmax = 0, CB_var_h = 0, SB_var_gmax = 0, SB_var_h = 0, PB_var_gmax = 0, PB_var_h = 0, variability = function(strain_parameter, var) { strain_parameter * 2^(seq(-var, var, length = length(strain_parameter))) } )add_strain_var( strain_parameter, CB_var_gmax = 0, CB_var_h = 0, SB_var_gmax = 0, SB_var_h = 0, PB_var_gmax = 0, PB_var_h = 0, variability = function(strain_parameter, var) { strain_parameter * 2^(seq(-var, var, length = length(strain_parameter))) } )
strain_parameter |
object of class |
CB_var_gmax |
the argument var for the function |
CB_var_h |
the argument var for the function |
SB_var_gmax |
the argument var for the function |
SB_var_h |
the argument var for the function |
PB_var_gmax |
the argument var for the function |
PB_var_h |
the argument var for the function |
variability |
function of which takes two arguments, i.e.
|
the value of strain_parameter with added strain variability
add_strain_var(new_strain_parameter(n_CB = 3, n_PB = 3, n_SB = 3), 2)add_strain_var(new_strain_parameter(n_CB = 3, n_PB = 3, n_SB = 3), 2)
a_0 potential added, and the possibility to simulate
multiple strains per functional groupThe rate equations, as published in the Bush et al 2017
doi:10.1093/clinchem/39.5.766 paper, but with forcing of
oxygen diffusivity a_0 potential added, and the possibility to simulate
multiple strains per functional group
bushplus_dynamic_model(t, state, parameters, log10a_forcing_func, ...)bushplus_dynamic_model(t, state, parameters, log10a_forcing_func, ...)
t |
The current time in the simulation |
state |
A vector containing the current (named) values of each state variable |
parameters |
An object of class |
log10a_forcing_func |
function to change oxygen diffusivity |
... |
not used. Needed to catch additional parameters. |
a list containing two elements, namely the rate of change of the
strains, and also the current values of oxygen diffusivity a.
This function contains events that can alter the state variables. In this event definition, densities are given a floor of 1 at every event.
event_definition_1( times, state, parms, log10a_forcing_func, noise_sigma, minimum_abundances )event_definition_1( times, state, parms, log10a_forcing_func, noise_sigma, minimum_abundances )
times |
the time points in the simulation when events happen |
state |
current state variable values |
parms |
object of class |
log10a_forcing_func |
function to change oxygen diffusivity |
noise_sigma |
value of variation added to |
minimum_abundances |
minimum abundances for CB, PB and SB. if the values
are lower, they will be set to these |
the state vector, i.e. a vector containing the state variables.
This function contains events that can alter the state variables. In this event definition, the abundances of all functional groups have 1 added to their density at each event.
event_definition_2( times, state, parms, log10a_forcing_func, noise_sigma, minimum_abundances )event_definition_2( times, state, parms, log10a_forcing_func, noise_sigma, minimum_abundances )
times |
the time points in the simulation when events happen |
state |
current state variable values |
parms |
object of class |
log10a_forcing_func |
function to change oxygen diffusivity |
noise_sigma |
value of variation added to |
minimum_abundances |
minimum abundances for CB, PB and SB. if the values
are lower, they will be set to these |
the state vector, i.e. a vector containing the state variables.
This function contains events that can alter the state variables. In this event definition, the abundances of all functional groups have 1 added to their density at each event.
event_definition_symmetric( times, state, parms, log10aO_forcing_func, log10aS_forcing_func, noise_sigma, minimum_abundances )event_definition_symmetric( times, state, parms, log10aO_forcing_func, log10aS_forcing_func, noise_sigma, minimum_abundances )
times |
the time points in the simulation when events happen |
state |
current state variable values |
parms |
object of class |
log10aO_forcing_func |
function to change oxygen diffusivity |
log10aS_forcing_func |
function. to change sulfide diffusivity |
noise_sigma |
value of variation added to |
minimum_abundances |
minimum abundances for CB, PB and SB. if the values
are lower, they will be set to these |
the state vector, i.e. a vector containing the state variables.
Can be used directly. Also can be used indirectly, via the function run_replication_ssfind_parameter()
get_final_states_a_N(x, parameter)get_final_states_a_N(x, parameter)
x |
A vector of oxygen diffusivity and initial conditions. |
parameter |
object of class |
A vector of steady states
Get the maximum of environmental conditions for which alternate stable states exist
get_hysteresis_max(up, down, a, threshold_diff)get_hysteresis_max(up, down, a, threshold_diff)
up |
State variable values as the environmental condition increases |
down |
State variable values as the environmental condition decreases |
a |
An environmental driver, here it is usually oxygen diffusivity |
threshold_diff |
Amount of different between up and down states to qualify as alternate stable states |
A numeric value, which is the extent of the range of conditions for which alternate stable states exist.
Get the minimum of environmental conditions for which alternate stable states exist
get_hysteresis_min(up, down, a, threshold_diff)get_hysteresis_min(up, down, a, threshold_diff)
up |
State variable values as the environmental condition increases |
down |
State variable values as the environmental condition decreases |
a |
An environmental driver, here it is usually oxygen diffusivity |
threshold_diff |
Amount of different between up and down states to qualify as alternate stable states |
A numeric value, which is the extent of the range of conditions for which alternate stable states exist.
Get the range of environmental conditions for which alternate stable states exist
get_hysteresis_range(up, down, a, threshold_diff)get_hysteresis_range(up, down, a, threshold_diff)
up |
State variable values as the environmental condition increases |
down |
State variable values as the environmental condition decreases |
a |
An environmental driver, here it is usually oxygen diffusivity |
threshold_diff |
Amount of different between up and down states to qualify as alternate stable states |
A numeric value, which is the extent of the range of conditions for which alternate stable states exist.
Get the total hysteresis of a system variable.
get_hysteresis_total(up, down)get_hysteresis_total(up, down)
up |
State variable values as the environmental condition increases |
down |
State variable values as the environmental condition decreases |
The total hysteresis.
Method described in Emancipator & Knoll (1993) Clinical Chemistry 39, 766-772, doi:10.1093/clinchem/39.5.766
get_nonlinearity(x, y)get_nonlinearity(x, y)
x |
Environmental driver state variable |
y |
Ecosystem state variable |
Measure of non-linearity
runsim_parameter from replication_ssfind_parameter
Extract the parameter to run an individual simulation from either a parameter set for a
steady state determination (run_replication_ssfind_parameter()) or the result set returned.
get_ssbyaN_parameter(result, i)get_ssbyaN_parameter(result, i)
result |
object of class |
i |
Index of the individual simulation to extract the parameter for |
an object of class replication_ssfind_parameter
The measures include non-linearity and hysteresis measures, of an ecosystem response to environmental change. Takes steady state data as the input.
get_stability_measures(ss_object, threshold_diff_log10scale, ...) ## S3 method for class 'replication_ssfind_result' get_stability_measures(ss_object, threshold_diff_log10scale = 3, ...) get_stability_measures_replication_ssfind_result(ss_object) ## S3 method for class 'temporal_ssfind_result' get_stability_measures(ss_object, threshold_diff_log10scale = 3, ...) get_stability_measures_temporal_ssfind_result(ss_object)get_stability_measures(ss_object, threshold_diff_log10scale, ...) ## S3 method for class 'replication_ssfind_result' get_stability_measures(ss_object, threshold_diff_log10scale = 3, ...) get_stability_measures_replication_ssfind_result(ss_object) ## S3 method for class 'temporal_ssfind_result' get_stability_measures(ss_object, threshold_diff_log10scale = 3, ...) get_stability_measures_temporal_ssfind_result(ss_object)
ss_object |
object from which to calculate the stability measures. At
the moment |
threshold_diff_log10scale |
Amount of different between up and down states to qualify as alternate stable states, on log10 scale |
... |
additional arguments for methods |
A data frame of stability measures of each state variable
Computes measures that make comparisons between collapse and recovery trajectory and between antagonistic environmental variables possible.
get_symmetry_measurements(res)get_symmetry_measurements(res)
res |
Results from |
returns a frame containing symmetry measures for each variable, such as hysteresis area, shift magnitudes or distance between TP, as well as TP themselfes with their according state.
Growth rate function on one resource X
growth1(X, g_max, k_X)growth1(X, g_max, k_X)
X |
Concentration of resource X |
g_max |
Maximum growth rate |
k_X |
Half saturation constant for resource X |
Growth rate
Growth rate function on two resources X and Y
growth2(X, Y, g_max, k_X, k_Y)growth2(X, Y, g_max, k_X, k_Y)
X |
Concentration of resource X |
Y |
Concentration of resource Y |
g_max |
Maximum growth rate |
k_X |
Half saturation constant for resource X |
k_Y |
Half saturation constant for resource Y |
Growth rate
Growth inhibition function
inhibition(X, h_X)inhibition(X, h_X)
X |
Concentration of substance X |
h_X |
Concentration of substance X at which the inhibition factor is 0.5 (i.e. the concerned rate is halved) |
Inhibition factor
Create CB strain parameter values
new_CB_strain_parameter(n = 1, values = "bush")new_CB_strain_parameter(n = 1, values = "bush")
n |
number of strains |
values |
Allowed values are:
|
object of class CB_strain_parameter. The object contains a
data.frame with 8 columns and n rows. The columns are:
columns:
strain_name: the name of the strain
g_max_CB:
k_CB_P:
h_SR_CB:
y_P_CB:
Pr_CB:
m_CB:
i_CB:
Create initial state of the system, selecting them from various preset options.
new_initial_state(n_CB = 1, n_PB = 1, n_SB = 1, values = "bush_anoxic_fig2ab")new_initial_state(n_CB = 1, n_PB = 1, n_SB = 1, values = "bush_anoxic_fig2ab")
n_CB |
number of CB strains |
n_PB |
number of PB strains |
n_SB |
number of SB strains |
values |
Allowed values are:
|
initial state vector aof the system s part of the strain_starter
Create PB strain parameter
new_PB_strain_parameter(n = 1, values = "bush")new_PB_strain_parameter(n = 1, values = "bush")
n |
number of strains |
values |
Allowed values are:
|
object of class PB_strain_parameter. The object contains a
data.frame with 8 columns and n rows. The columns are:
columns:
strain_name: the name of the strain
g_max_PB:
k_PB_SR:
k_PB_P:
h_O_PB:
y_SR_PB:
y_P_PB:
m_PB:
i_CB:
run_replication_ssfind_parameter() to run such a set of simulations.Create parameter set to run a set of simulations to find stable states. Is passed to the function run_replication_ssfind_parameter() to run such a set of simulations.
new_replication_ssfind_parameter(...)new_replication_ssfind_parameter(...)
... |
named parameter for the simulation to be set. An error will be raised, if they are not part of the parameter set. |
object of the class replication_ssfind_parameter. This object of class
replication_ssfind_parameter is identical to an object of class runsim_parameter
plus an additional field ss_expt which contains a data.frame with
columns named
N_CB
N_PB
N_SB
a_0
which contain the initial states. If more than one strain is present, the same
initial state is assumed for each.
replication_ssfind_result which is returned by the
function run_replication_ssfind().Create object of type replication_ssfind_result which is returned by the
function run_replication_ssfind().
new_replication_ssfind_results(parameter, result)new_replication_ssfind_results(parameter, result)
parameter |
object of class |
result |
a |
object of the class replication_ssfind_result. This object of class
replication_ssfind_result is identical to an object of class replication_ssfind_parameter
plus an additional field result which contains the result of the simulation.
Create parameter set with which to run a simulation.
new_runsim_parameter(...)new_runsim_parameter(...)
... |
named parameter for the simulation to be set. An error will be raised, if they are not part of the parameter set. |
parameter object of the class runsim_parameter.
The object contains the following elements:
dynamic_model : the dynamic model to be used. At the moment, only bushplus_dynamic_model is implemented. Fur further info, see the documen tation of bushplus_dynamic_model.
event_definition : A function which alters the state variables. At the moment only event_definition_1() is included. User defined functions with the same signature can be used.
strain_parameter : object of class strain_parameter as returned by new_strain_parameter()
event_interval : interval, in timesteps, in which the event occurs"
noise_sigma : value of variation added to SO, SR, O, and P during the event
minimum_abundances : Minimum abundances. Smaller abundances will be set to this value duting event_definition_1().
sim_duration : duration of the simulation
sim_sample_interval: interval, at which the simulation will be sampled
log10a_series : A vector of values of log10 oxygen diffusivity parameter at which stable states will be found.
asym_factor: : For symmetric simulations only: enables manipulating aS forcing in asymmetric manner to decrease (<1) or increase (>1) stress on cyanobacteria.
solver_method : Used for the solver. Default is "radau". For other options, see the documentatioom of odeSolve::ode.
runsim_result which is returned by the
function run_sim().Create object of type runsim_result which is returned by the
function run_sim().
new_runsim_results(parameter, result)new_runsim_results(parameter, result)
parameter |
object of class |
result |
a |
object of the class runsim_result. This object of class
runsim_result is identical to an object of class runsim_parameter
plus an additional field result which contains the result of the simulation.
SB_strain_parameter objectCreate new SB_strain_parameter object
new_SB_strain_parameter(n = 1, values = "bush")new_SB_strain_parameter(n = 1, values = "bush")
n |
number of strains |
values |
Allowed values are:
|
object of class SB_strain_parameter. The object contains a
data.frame with 8 columns and n rows. The columns are:
columns:
strain_name: the name of the strain
g_max_SB:
k_PB_SO:
k_SB_P:
h_O_SB:
y_SO_SB:
y_PB_SB:
m_SB:
i_SB:
strain_parameter
Creates a set of parameters and starting conditions for a simulation. This function assists with the initialization of a simulation, by providing various reference sets of parameter values and initial conditions. The default is to create the parameter set used in Bush et al. (2017) doi:10.1093/clinchem/39.5.766 and the anoxic favorable initial conditions used in the simulations for Figure 2a&b of that publication.
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" )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" )
n_CB |
number of CB strains, default 1. |
values_CB |
Allowed values are:
|
n_PB |
number of PB strains, default 1. |
values_PB |
Allowed values are:
|
n_SB |
number of SB strains, default 1. |
values_SB |
Allowed values are:
|
values_other |
Allowed values are:
|
values_initial_state |
values to be used for initial values or
|
Object of class strain_parameter. The object contains the following fields:
The additional parameters are:
CB: strain parameter from Cyano Bacteria
PB: strain parameter for the Phototrophic bacteria
SB: strain parameter for the Sulphur bacteria
a_S: substrate diffusivity of sulphur <- 0.001
a_O: substrate diffusivity of oxygen <- 8e-4
a_P: substrate diffusivity of phosphorous <- 0.01
back_SR: background substrate concentration of XXX <- 300
back_SO: background substrate concentration of XXX <- 300
back_O: background substrate concentration of oxygen <- 300
back_P: background substrate concentration of phosphorous <- 9.5
c: <- 4e-5
temporal_ssfind_result which is returned by the
function run_temporal_ssfind_parameter().Create object of type temporal_ssfind_result which is returned by the
function run_temporal_ssfind_parameter().
new_temporal_ssfind_results(parameter = NULL, result)new_temporal_ssfind_results(parameter = NULL, result)
parameter |
not used - for consistency. |
result |
a |
the result object.
This is a convenience function to plot the dynamics of a model run, with strains within functional groups displayed.
plot_dynamics(simulation_result, every_n = 1)plot_dynamics(simulation_result, every_n = 1)
simulation_result |
Object returned by the run_simulation function |
every_n |
Plot data of every other n sample. |
returns the ggplot object of the plot. If it is assigned to a variable, the plot needs to be plotted, otherwise it is plotted.
This is a convenience function to plot the dynamics of a model run particularly for a symmetric model with strains within functional groups displayed.
plot_dynamics_symmetric(simulation_result, every_n = 1, plot_a = FALSE)plot_dynamics_symmetric(simulation_result, every_n = 1, plot_a = FALSE)
simulation_result |
Object returned by the run_simulation function |
every_n |
Plot data of every other n sample. |
plot_a |
if |
returns the ggplot object of the plot. If it is assigned to a variable, the plot needs to be plotted, otherwise it is plotted.
get_symmetry_measures()
Visualizes the symmetry measures obtained by the function get_symmetry_measures()
plot_symmetry_measures(res, species)plot_symmetry_measures(res, species)
res |
Results from |
species |
Environmental varibale of which the symmetry measures shoudl be plotted |
ggplot object, with shifts visualized as arrows and hysteresis area ribbon
This is a convenience function to plot the dynamics of a model run particularly for a symmetric model with strains within functional groups displayed.
plot_temporal_ss(temporal_results)plot_temporal_ss(temporal_results)
temporal_results |
Results from |
returns the ggplot object of the plot. If it is assigned to a variable, the plot needs to be plotted, otherwise it is plotted.
get_symmetry_measurements()
Visualizes the symmetry of antagonistic collapse or recovery trajectories, symmetry measures of shift
are computed by the function get_symmetry_measurements()
plot_trajectory_symmetry( res, trajectory = "recovery", typ = "substrate", plot_log10 = FALSE )plot_trajectory_symmetry( res, trajectory = "recovery", typ = "substrate", plot_log10 = FALSE )
res |
Results from |
trajectory |
Trajectory to be plotted. Either 'collapse' or 'recovery'.'. |
typ |
Type of environmental variable. Either 'substrate' or 'bacteria'. |
plot_log10 |
Logical. If |
ggplot object, with antagonistic shifts of related trajectory visualized as arrows.
get_symmetry_measurements()
Visualizes the symmetry of upards and downwards trajectories in the same figure,
symmetry measures of shift are computed by the function get_symmetry_measurements()
plot_trajectory_symmetry_compact(res, typ = "substrate", plot_log10 = FALSE)plot_trajectory_symmetry_compact(res, typ = "substrate", plot_log10 = FALSE)
res |
Results from |
typ |
Type of environmental variable. Either 'substrate' or 'bacteria'. |
plot_log10 |
Logical. If |
ggplot object, with antagonistic shifts of related trajectory visualized as arrows.
Function to get the steady states for combinations of a (oxygen diffusivity) and initial states.
This function is multithreaded, and the value of mc.cores determines the umber of parallel threads.
run_replication_ssfind(parameter, mc.cores = getOption("mc.cores", 0))run_replication_ssfind(parameter, mc.cores = getOption("mc.cores", 0))
parameter |
object of class |
mc.cores |
the number of cores to be used. If 0, the old sequential version is used.
The default is read from the option |
Processed data about steady states
This function takes the parameter object and runs a simulation based on these.
It returns an object of class runsim_result which contains an additional
entry, i.e. result which contains the results of the simulation. The
simulation can be re-run using the returned object as input parameter.
run_simulation(parameter)run_simulation(parameter)
parameter |
an object of class |
an object of class runsim_result, obtained from running the
simulation as defined in parameter.'
This function takes the parameter object and runs a simulation based on these.
It returns an object of class runsim_result which contains an additional
entry, i.e. result which contains the results of the simulation. The
simulation can be re-run using the returned object as input parameter.
run_simulation_symmetric(parameter)run_simulation_symmetric(parameter)
parameter |
an object of class |
an object of class runsim_result, obtained from running the
simulation as defined in parameter.'
parameter$sim_duration
while keeping the length of parameter$log10_series doesn't change response dynanimcs,
stable states have been found.Used to find the stable states for a parameter set by increasing and decreasing
the oxygen diffusivity in a stepwise fashion. If increasing parameter$sim_duration
while keeping the length of parameter$log10_series doesn't change response dynanimcs,
stable states have been found.
run_temporal_ssfind(parameter)run_temporal_ssfind(parameter)
parameter |
an object of class |
A data frame of final states and oxygen diffusivity values
Calls the function run_temporal_ssfind for each parameter set.
run_temporal_ssfind_experiment( parameter, var_expt, total_initial_abundances, cores = 1 )run_temporal_ssfind_experiment( parameter, var_expt, total_initial_abundances, cores = 1 )
parameter |
an object of class |
var_expt |
An object that describes different levels of diversity that should be examined. This object is not created in the 'microxanox' package. |
total_initial_abundances |
An object containing the total abundance in each functional group. |
cores |
Number of cores to use. If more than one, then |
an tibble object containing the stable state result, as well as the simulation parameters.
parameter$sim_duration
while keeping the length of parameter$log10_series doesn't change response dynanimcs,
stable states have been found.Used to find the stable states for a parameter set by increasing and decreasing
the oxygen diffusivity in a stepwise fashion. If increasing parameter$sim_duration
while keeping the length of parameter$log10_series doesn't change response dynanimcs,
stable states have been found.
run_temporal_ssfind_symmetric(parameter)run_temporal_ssfind_symmetric(parameter)
parameter |
an object of class |
A data frame of final states, oxygen and sulfide diffusivity values
The set_diffusivities enables manipulating aS forcing in asymmetric manner to decrease (<1) or increase (>1) stress on cyanobacteria. If the vector in
parameter$log10a_series has multiple local minima, maxima respectively, it is directly mirrored. ATTENTION: Mean value of the series has to be the axis at which the
vector is mirrored!
set_diffusivities(param, asym_factor = 1)set_diffusivities(param, asym_factor = 1)
param |
:parameter set of type |
asym_factor |
:asym_factor to update |
updated parameterset
run_temporal_ssfind_experiment
Should only be used internally.Convenience function for setting initial states during the process of finding
stable states using the function run_temporal_ssfind_experiment
Should only be used internally.
set_temporal_ssfind_initial_state( p, initial_total_CB, initial_total_PB, initial_total_SB )set_temporal_ssfind_initial_state( p, initial_total_CB, initial_total_PB, initial_total_SB )
p |
parameters for the simulation |
initial_total_CB |
Total initial abundance of CB |
initial_total_PB |
Total initial abundance of PB |
initial_total_SB |
Total initial abundance of SB |
The passed parameter object but with initial states set (overwritten by the new ones)
The rate equations, modified from publication of Bush et al 2017 doi:10.1093/clinchem/39.5.766. Equations of antagonistic environmental variables are identical, while the oxygen and sulfide diffusivities are forced in antisymmetric manner. Possibility to simulate multiple strain per functional group as published in doi:10.1111/ele.14217 by Limberger et al 2023 remains.
symmetric_bushplus_dynamic_model( t, state, parameters, log10aO_forcing_func, log10aS_forcing_func, ... )symmetric_bushplus_dynamic_model( t, state, parameters, log10aO_forcing_func, log10aS_forcing_func, ... )
t |
The current time in the simulation |
state |
A vector containing the current (named) values of each state variable |
parameters |
An object of class |
log10aO_forcing_func |
function to change oxygen diffusivity |
log10aS_forcing_func |
function to change sulfide diffusivity |
... |
not used. Needed to catch additional parameters. |
a list containing two elements, namely the rate of change of the
strains, and also the current values of oxygen diffusivity aO, as well as sulfide diffusivity aS.