vignettes/rescomp.Rmd
rescomp.Rmd
rescomp is an R package for generating, simulating and visualising ODE models of consumer-resource interactions. In essence, it is a consumer-resource modelling focused interface to the canonical deSolve package.
All rescomp
models take the general form,
\[ \frac{dN_{i}}{dt} = N_{i}(\mu_{i}(R_{1},R_{2},...)-m) \,, \]
\[ \frac{dR_{j}}{dt} = \Psi_{j} (R_{j}) - \sum_{i = 1}^{n}Q_{ij} \mu_{ij}(R_j) N_i \,, \]
where \(N_{i}\) is the population density of consumer \(i\), \(R_{j}\) is the density/concentration of resource \(j\), \(\mu_{i}()\) is the per capita consumer functional response of consumer \(i\), \(m\) is the per capita mortality rate (constant or intermittent), \(\Psi_{j}(R_j)\) is the resource supply function, and \(Q_{ij}\) is the resource quota of consumer \(i\) on resource \(j\) (the amount of resource per unit consumer).
From this general form, different model formulations are distinguished by: i) the number of consumers/resources; ii) the form of the consumer functional response; iii) the mode of resource supply; iv) the type of resource; and v) any non-autonomous behaviour including time dependent model parameters and/or instantaneous changes in consumer/resource density (e.g. in batch transfer).
The consumer growth function can take one of three forms:
\[ \mu_{ij}(R_{j}) = a_{ij}R_{j} \,, \]
where \(a_{ij}\) is a resource consumption constant.
\[ \mu_{ij}(R_{j}) = \mu_{max_{ij}}\frac{R_{j}}{k_{ij} + R_{j}} \,, \]
where \(\mu_{max_{ij}}\) is the maximum growth rate and \(k_{ij}\) is the half saturation constant for consumer \(i\) on resource \(j\).
\[ \mu_{ij}(R_{j}) = \mu _{max_{ij}}\frac{R_{j}^2}{k_{ij}^2 + R_{j}^2} \,. \]
An alternative parameterisation to all three functional responses is to incorporate a consumption efficiency term, \(c_{ij}\). The type I model becomes \(\mu_{ij}(R_{j}) = a_{ij}c_{ij}R_{j}\), and for type II and III, \(a_{ij}c_{ij}\) is substituted for \(\mu _{max_{ij}}\). NB. this functional form is mutually exclusive with the specification of resources quotas in the resource equations, i.e. the \(Q_{ij}\)s are dropped from the resource equations.
The resource supply function, \(\Psi_{j}(R_j)\), can also take several forms. Specifically, either the resources are biological and grow logistically,
\[ \Psi_{j} (R_{j}) = r_{j}R_{j}\left(1-\frac{R_{j}}{K_{j}}\right) \,, \]
where \(r_j\) is the resource intrinsic rate of increase and \(K_j\) is the resource carrying capacity;
Or the resources are supplied to the systems at a fixed concentration and rate (as in a chemostat),
\[ \Psi_{j} (R_{j}) = d(S_{j} - R_{j}) \,, \]
where \(d\) represents the flux of resources into and out of the system;
Or the resources are pulsed intermittently (e.g., in serial transfer - see examples below).
In the case of multiple resources, each resource is either treated as essential to consumer growth following Leibig’s law of the minimum, in which case,
\[ \mu_{i}(R_{1}, R_{2},...,R_{n}) = min(\mu_{i}(R_{1}), \mu_{i}(R_{2}),..., \mu_{i}(R_{n})) \,; \]
or substitutable such that:
\[ \mu_{i}(R_{1}, R_{2},...,R_{n}) = \mu_{i}(R_{1}) + \mu_{i}(R_{2}) + ... + \mu_{i}(R_{n}) \,. \]
rescomp
The main user function in rescomp is spec_rescomp()
, which facilitates the definition and parameterisation of a consumer-resource model and the specification of simulation parameters. Function arguments include (but are not limited to):
See ?spec_rescomp
for all argument options.
Examples demonstrated in this vignette include:
One consumer (type 2) and one logistically growing resource
One consumer (type 2) and one continuously supplied resource (i.e. chemostat)
Two consumers (type 1 and type 2) and one logistically growing resource
Two consumers (type 1) and two substitutable resources in a chemostat
Two consumers (type 1) and two essential resources in a chemostat
Two consumers (type 2) and one externally pulsed resource
Two consumers (type 2) with time dependent consumption parameters and one continuously supplied resource
Three consumers (type 1) on three essential resources generating an intransitive loop (i.e. rock paper scissors dynamics).
Additional examples illustrating other options and functionality can be found in the vignette “Reproducing studies in resource competition with rescomp”, including delayed consumer introduction times and specifying a model with consumer resource efficiencies rather than resource quotas.
pars <- spec_rescomp(
spnum = 1,
resnum = 1,
funcresp = "type2",
mumatrix = matrix(0.12),
chemo = FALSE,
resspeed = 3,
resconc = 2,
totaltime = 500
)
#> Model properties
#> * 1 consumer(s) and 1 resource(s)
#> * Consumers have type 2 functional responses
#> * Resources grow logistically
#> * Mortality is continuous
#>
#> Simulation properties
#> * Simulation time: 500 time steps
#> * Init state: consumer(s) = [10], resource(s) = [2]
Visualise functional responses with plot_funcresp()
. Faceted by resources and time dependent parameterisations (when relevant).
plot_funcresp(pars, maxx = 2)
The model can then be simulated with sim_rescomp()
. sim_rescomp()
is a wrapper to deSolve::ode()
that takes the output from spec_rescomp()
as its main argument.
m1 <- sim_rescomp(pars)
plot_rescomp()
produces a plot of the output dynamics (both consumers and resources).
plot_rescomp(m1)
Note the total simulation time and initial state variables specified with spec_rescomp()
can be overidden in sim_recomp()
with the arguments times
and y
, respectively. This will nevertheless print a message to the console to check it is intentional. See ?time_vals()
for convenient specification of simulation time (and pulsing when used).
m1 <- sim_rescomp(pars, y = c(30000, 1), times = time_vals(200))
#> totaltime in pars will be overidden
#> cinit in pars will be overidden
plot_rescomp(m1)
pars <- spec_rescomp(
spnum = 1,
resnum = 1,
funcresp = "type2",
mumatrix = matrix(0.12),
resspeed = 0.03,
resconc = 2,
totaltime = 300
)
#> Model properties
#> * 1 consumer(s) and 1 resource(s)
#> * Consumers have type 2 functional responses
#> * Resource supply is continuous (e.g. chemostat)
#> * Mortality is continuous
#>
#> Simulation properties
#> * Simulation time: 300 time steps
#> * Init state: consumer(s) = [10], resource(s) = [2]
m1 <- sim_rescomp(pars)
plot_rescomp(m1)
pars <- spec_rescomp(
spnum = 2,
resnum = 1,
funcresp = c("type1", "type2"),
mumatrix = matrix(c(0.35,0.05),
nrow = 2,
ncol = 1,
byrow = TRUE),
kmatrix = matrix(c(1, 0.015),
nrow = 2,
ncol = 1,
byrow = TRUE),
chemo = FALSE,
resspeed = 1,
resconc = 0.2,
totaltime = 2000
)
#> Model properties
#> * 2 consumer(s) and 1 resource(s)
#> * Consumers have type 1 or type 2 functional responses
#> * Resources grow logistically
#> * Mortality is continuous
#>
#> Simulation properties
#> * Simulation time: 2000 time steps
#> * Init state: consumer(s) = [10, 10], resource(s) = [0.2]
plot_funcresp(pars, maxx = 0.2)
m1 <- sim_rescomp(pars)
plot_rescomp(m1)
Owing to the trade-off in functional responses (gleaner-opportunist) and resource fluctuations generated by N2, the two consumers are able to coexist.
pars <- spec_rescomp(
spnum = 2,
resnum = 2,
funcresp = "type1",
mumatrix = matrix(c(0.09,0.04,
0.05,0.08),
nrow = 2,
ncol = 2,
byrow = TRUE),
resspeed = 0.03,
resconc = 1,
mort = 0.03,
essential = FALSE,
totaltime = 1000
)
#> Model properties
#> * 2 consumer(s) and 2 resource(s)
#> * Consumers have type 1 functional responses
#> * Resources are substitutable
#> * Resource supply is continuous (e.g. chemostat)
#> * Mortality is continuous
#>
#> Simulation properties
#> * Simulation time: 1000 time steps
#> * Init state: consumer(s) = [10, 10], resource(s) = [1, 1]
plot_funcresp(pars, maxx = 1)
m1 <- sim_rescomp(pars)
plot_rescomp(m1)
Resource partitioning, where each consumer is a better competitor (lower R*) for a different resource, results in coexistence.
Using pars$essential = TRUE
switches the previous parameteristion to essential resources while keeping all else equal.
pars$essential = TRUE
m1 <- sim_rescomp(pars)
plot_rescomp(m1)
There is no coexistence now because each species needs to have a larger impact on the resource that is most limiting its growth. For ‘consumer 1’ it is ‘resource b’, whereas for ‘consumer 2’ it is ‘resource a’. We can adjust consumer resource impacts via the resource quota matrix.
pars <- spec_rescomp(
spnum = 2,
resnum = 2,
funcresp = "type1",
mumatrix = matrix(c(0.09, 0.04,
0.05, 0.08),
nrow = 2,
ncol = 2,
byrow = TRUE),
qmatrix = matrix(c(0.001, 0.005,
0.005, 0.001),
nrow = 2,
ncol = 2,
byrow = TRUE),
resspeed = 0.03,
resconc = 1,
mort = 0.03,
essential = TRUE
)
#> Model properties
#> * 2 consumer(s) and 2 resource(s)
#> * Consumers have type 1 functional responses
#> * Resources are essential
#> * Resource supply is continuous (e.g. chemostat)
#> * Mortality is continuous
#>
#> Simulation properties
#> * Simulation time: 1000 time steps
#> * Init state: consumer(s) = [10, 10], resource(s) = [1, 1]
m1 <- sim_rescomp(pars)
plot_rescomp(m1)
External resource pulsing is handled using an auxiliary event function which is generated by the call to spec_rescomp()
and passed to sim_rescomp()
. The size of the resource pulse is specified as a numeric with the argument respulse
in the call to spec_rescomp()
. To prevent any other mode of resource supply, resspeed
should be set to zero. The frequency of pulsing is set with the pulsefreq
argument.
pars <- spec_rescomp(
spnum = 2,
resnum = 1,
funcresp = "type2",
mumatrix = matrix(c(0.7,0.05),
nrow = 2,
ncol = 1,
byrow = TRUE),
kmatrix = matrix(c(2, 0.015),
nrow = 2,
ncol = 1,
byrow = TRUE),
resspeed = 0, # set to zero for no additional resource supply
resconc = 0.2,
respulse = 0.3,
pulsefreq = 100 # resource pulse size
)
#> Model properties
#> * 2 consumer(s) and 1 resource(s)
#> * Consumers have type 2 functional responses
#> * Resource supply is pulsed only
#> * Mortality is continuous
#>
#> Simulation properties
#> * Simulation time: 1000 time steps
#> * Resources pulsing every 100 timesteps
#> * Init state: consumer(s) = [10, 10], resource(s) = [0.2]
plot_funcresp(pars, maxx = max(pars$respulse))
m1 <- sim_rescomp(pars)
plot_rescomp(m1)
Serial transfer in batch culture differs from the above scenario in that resource pulses are typically accompanied by large instantaneous changes in population density. For example, in a serial transfer experiment with bacteria, an experimenter might sample a fraction (e.g. 10%) of the community every 24 hours and inoculate that fraction intro fresh media. As such, the consumer population density immediately following a transfer event will be equal to: \((1-M)N_{i}\), where \(M\) is instantaneous mortality fraction given by mortpulse
.
Simulating serial transfer in batch culture is implemented by:
batchtrans = TRUE
in spec_rescomp()
.mortpulse
. This is the fraction of the consumer(s) population instantaneous ‘killed’ at set intervals. Note the transfer fraction = 1 - mortpulse
.pulsefreq
argument (coincident with resource pulse intervals).Note that the effect of setting batchtrans = TRUE
is to enforce fractional sampling of the resource/media in addition to the consumers. As such, immediately after a transfer event the resource concentration will be equal to: \((1-M)R_{j} + MS_{j}\), where \(R\) is the resource concentration prior to transfer, and \(S\) is the resource pulse size given by respulse
. To inhibit this behaviour, i.e. to stop fractional sampling of the resource but maintain intermittent mortality use batchtrans = FALSE
but still assign a non-zero value for mortpulse
.
pars <- spec_rescomp(
spnum = 2,
resnum = 1,
funcresp = "type2",
mumatrix = matrix(c(0.2,0.2),
nrow = 2,
ncol = 1,
byrow = TRUE),
kmatrix = matrix(c(0.3, 0.2),
nrow = 2,
ncol = 1,
byrow = TRUE),
resspeed = 0,
resconc = 0.6,
respulse = 1,
mort = 0, # set to zero to make the transfer fraction the effective mortality rate
mortpulse = 0.8,
batchtrans = TRUE,
pulsefreq = 100,
)
#> Model properties
#> * 2 consumer(s) and 1 resource(s)
#> * Consumers have type 2 functional responses
#> * Resource supply is pulsed only
#> * Mortality intermittent
#>
#> Simulation properties
#> * Simulation time: 1000 time steps
#> * Resources pulsing and intermittent mortality every 100 timesteps
#> * Init state: consumer(s) = [10, 10], resource(s) = [0.6]
plot_funcresp(pars, maxx = max(pars$respulse))
m1 <- sim_rescomp(pars)
plot_rescomp(m1)
Note, to reproduce the dynamics of a single batch culture experiment (i.e. without serial transfer) as potentially measured via optical density in a plate reader (negligible or unobserved mortality), simply set both resspeed
and mort
to zero.
pars <- spec_rescomp(resspeed = 0, mort = 0, totaltime = 500)
#> Model properties
#> * 1 consumer(s) and 1 resource(s)
#> * Consumers have type 1 functional responses
#> * Resources are not supplied?!
#> * No mortality
#>
#> Simulation properties
#> * Simulation time: 500 time steps
#> * Init state: consumer(s) = [10], resource(s) = [1]
m1 <- sim_rescomp(pars)
plot_rescomp(m1)
Using time dependent consumption parameters is appropriate for capturing basic consumer-environment interactions (e.g. differential responses to temperature, pH, antibiotics etc.). Time varying parameters can be implemented in with forcing functions. rescomp provides a convenient means of implementing forcing functions in the background. Two additional arguments are needed in calls to spec_rescomp()
: timepars
and timeparfreq
. timepars
should be set to TRUE
; timeparfreq
adjusts the frequency of switching between ‘environmental states’.
In addition, if time_pars = TRUE
, the mumatrix argument in spec_rescomp
expects a list containing two matrices, reflecting two different environmental states. At this stage rescomp
only handles two states with three interpolation options: i) equal time allotted to each state (default); ii) linear (tpinterp = "lin"
) interpolation between the two states; iii) sinusoidal (tpinterp = "sine"
) interpolation between states.
pars <- spec_rescomp(
spnum = 2,
resnum = 1,
funcresp = "type2",
mumatrix = list(matrix(c(0.3,0.1),
nrow = 2,
ncol = 1,
byrow = TRUE),
matrix(c(0.2,0.5),
nrow = 2,
ncol = 1,
byrow = TRUE)),
kmatrix = matrix(c(0.1, 0.1),
nrow = 2,
ncol = 1,
byrow = TRUE),
chemo = TRUE,
resspeed = 0.03,
resconc = 0.1,
timepars = TRUE,
timeparfreq = 50,
totaltime = 2000
)
#> Model properties
#> * 2 consumer(s) and 1 resource(s)
#> * Consumers have type 2 functional responses
#> * Resource supply is continuous (e.g. chemostat)
#> * Mortality is continuous
#> * Time dependent parameters with instantaneous switching every 50 timesteps
#>
#> Simulation properties
#> * Simulation time: 2000 time steps
#> * Init state: consumer(s) = [10, 10], resource(s) = [0.1]
plot_funcresp(pars, maxx = max(pars$resconc))
m1 <- sim_rescomp(pars)
plot_rescomp(m1)
The two consumers coexist due species-specific environmental responses leading to a temporal storage effect.
pars <- spec_rescomp(
spnum = 3,
resnum = 3,
funcresp = "type1",
mumatrix = matrix(c(0.11, 0.07, 0.045,
0.05, 0.11, 0.07,
0.07, 0.047, 0.1),
nrow = 3,
ncol = 3,
byrow = TRUE),
qmatrix = matrix(c(0.2, 0.7, 0.4,
0.4, 0.2, 0.7,
0.7, 0.4, 0.2),
nrow = 3,
ncol = 3,
byrow = TRUE),
essential = TRUE,
chemo = FALSE,
mort = 0.015,
resspeed = 10,
resconc = 0.5,
totaltime = 2000
)
#> Model properties
#> * 3 consumer(s) and 3 resource(s)
#> * Consumers have type 1 functional responses
#> * Resources are essential
#> * Resources grow logistically
#> * Mortality is continuous
#>
#> Simulation properties
#> * Simulation time: 2000 time steps
#> * Init state: consumer(s) = [10, 10, 10], resource(s) = [0.5, 0.5, 0.5]
plot_funcresp(pars, maxx = max(pars$resconc))
m1 <- sim_rescomp(pars)
plot_rescomp(m1)