In flee-group/flume: A Model For Fluvial Meta-Ecosystems

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  error = TRUE
)

knitr::opts_chunk$set(fig.width = 5.5, fig.height = 5.5, collapse = TRUE, comment = "##")

library(flume)

Here we demonstrate the model using a simple simulation with generated data. There are several essential components:

1. A metacommunity object, which describes all the possible species that can occur in the simulation, gives their niches, and includes the pairwise competition strength.
2. A river_network, which describes the topology of the network and gives other important attributes, such as discharge.
3. State variables; these are (mostly) dynamic attributes describing the state of the model at a given time.
4. Boundary conditions: input, both to the biological community (via immigration from the surrounding landscape) and to the resources (via terrestrial runoff, groundwater intrusion, etc).

1. The metacommunity

The metacommunity is an object created by flume that describes all of the species that can possibly occur in the model. It also describes the possible resources and other habitat variables, their range of possible values, how the resources behave, the species' niches, how species compete for resources, and species' dispersal abilities.

Here we generate a metacommunity using a scenario. It is also possible to have finer control over how the metacommunity is modelled; see the "metacommunities" vignette for more details.

niche_args = list(scale_c = 2e-6, scale_e = 2e-7)
dispersal_args = list()
mc = metacommunity(nsp = 3, nr = 1, niches = niches_uniform, niche_args = niche_args, 
  dispersal_args = dispersal_args)
plot(mc)

The niches_uniform scenario spreads the resources out along a gradient, with species' optima uniformly spread along this gradient. We set the scale of the colonisation and extinction function (here in units of reaches/second), and leave the dispersal parameters blank for now. We can see the niche parameters used with the niche_par function:

cbind(
  niche_par(mc, "location"),
  niche_par(mc, "sd")
)

2. The river network

The main component of the river network is a weighted digraph; i.e., the topological relationships between locations. The network graph in flume uses equal-length stream reaches as nodes, while the edges describe which reaches are neighbours of each other. The simplest way to define a network is using an attribute list and a weighted adjacency matrix.

Discharge and cross-sectional area are required attributes (although area can be estimated from Q if omitted). Here, Q grows by 1 each time we move downstream (i.e., when we traverse an edge, representing lateral input), and we add reaches together at confluences. We define node 1 to be the most downstream node (this is purely an arbitrary choice). We define a river network with two equal subcatchments and a mainstem.

Q1 = Q2 = c(13,12,8,7,1,4,1,1,2,1)
Q3 = c(31,30,29,28,27)

The adjacency matrix describes the spatial relationships among nodes. The matrix must be square, with one row and column per node. For two sites i and j, adj[i,j] will be nonzero if i is upstream of j. The layout matrix shows where to plot each vertex; this is optional but can be useful for nicer plots.

It should be clear that creating a non-trivial network by hand is cumbersome! Fortunately there is a way to do this for a real river using a digital elevation model; see the network_import vignette for more details.

## matrix for the two subcatchments
nsites = length(Q1)
adj = matrix(0, nrow=nsites, ncol=nsites)
adj[2,1] = adj[3,2] = adj[9,2] = adj[4,3] = adj[5,4] = adj[6,4] = adj[7,6] = adj[8,6] = 
  adj[10,9] = 1
adj2 = adj

# layout is optional, and controls where to place the nodes
# here we just use values to make the network look somewhat nice
# For a real river, x-y spatial coordinates for each reach chould be used
layout1 = matrix(c(0,0,0,1,-0.5,2,-0.5,3,-1,4,0,4,-0.5,5,0.5,5,0.5,2,1,3), byrow=TRUE, nrow=nsites)
layout2 = layout1
layout2[,1] = layout2[,1] + 3

rownames(adj) = colnames(adj) = paste("c1", 1:nrow(adj), sep="_")
rownames(adj2) = colnames(adj2) = paste("c2", 1:nrow(adj2), sep="_")
nsites3 = length(Q3)
adj3 = matrix(0, nrow=nsites3, ncol=nsites3)
adj3[2,1] = adj3[3,2] = adj3[4,3] = adj3[5,4] = 1
rownames(adj3) = colnames(adj3) = paste("c3", 1:nrow(adj3), sep="_")
layout3 = matrix(c(rep(1.5, nsites3), (-nsites3):(-1)), ncol=2)

# now create an adjacency matrix for all subcatchments together
adj_all = matrix(0, nrow = nsites+nsites+nsites3, ncol = nsites+nsites+nsites3)
rownames(adj_all) = colnames(adj_all) = c(rownames(adj), rownames(adj2), rownames(adj3))

# add the links from the other matrices
adj_all[rownames(adj), colnames(adj)] = adj
adj_all[rownames(adj2), colnames(adj2)] = adj2
adj_all[rownames(adj3), colnames(adj3)] = adj3

# link catchments 1 and 2 to catchment 3
adj_all['c1_1', 'c3_5'] = adj_all['c2_1', 'c3_5'] = 1

network = river_network(adjacency = adj_all, discharge = c(Q1, Q2, Q3), 
  layout = rbind(layout1, layout2, layout3))

By default, plotting the network will show the edges weighted by discharge. Node ('vertex') colours can be set by vertex.color; I use a simple gradient following the order of sites, because this is also how I will set up my resource gradients (in 3. Setting up the model below). Edge weight is proportional to the discharge between sites.

plot(network)

3. Starting conditions

We will now set the starting resource conditions and species distributions. For this example, we have 3 species along a resource gradient (where the resource values are between 0 and 1). We set up the network here so that there is a spatial gradient in resources, with high concentration in one subcatchment, low in the other, and intermediate values in the mainstem. We then use an equilibrium scenario to establish initial species distributions.

st0 = rep(0.4, nrow(adj_all))

## headwaters at opposite sides of the gradient, then gradually come together moving downstream
st0[c(5,7,8,10)] = 0.85
st0[c(5,7,8,10)+10] = 0.1
st0[c(3,4,6,9)] = 0.75
st0[c(3,4,6,9)+10] = 0.15
st0[c(1:2)] = 0.7
st0[c(1:2)+10] = 0.2
st0 = matrix(st0, ncol=1)

state(network, "resources") = st0
vlc = c(rep("white", 10), rep("black", 10), rep("#444444", 5))
plot(network, vertex.label.color = vlc)

sp0 = community_equilibrium(network, mc)
state(network, "species") = sp0
plot(network, "species")

3. Setting up the model

The model is initialised with the flume() function. Here, we should also specify initial values for the two state variables (site by resource concentration, and site by species).

The state variables are a matrix, one row per node in the river network, with each column representing either resources or species. The resource state variable can be accessed with state(network), and the species with site_by_species(network).

By default, the boundary condition for resources will use identical values as the initial state, and the boundary condition for species is set to zero. These can be changed using the stb and spb arguments.

We use a scenario for the initial site by species matrix that places species at any site where they would be expected to be found at equilibrium in the absence of competition. This requires the initial state of the network to be specified in advance. Scenarios for resources are coming in a future version; for now we simply set an upstream to downstream gradient.

model = flume(mc, network, sp0, st0)

4. Running the model

Now we are prepared to run the simulation. Here we run for 400 time steps and then explore some of the output options. In this case, all species survive the entire simulation, and occupancy seems to fluctuate around the same(ish) values for all species.

set.seed(123)
model_run = run_simulation(model, nt = 400)

## occupancy plot
plot(model_run)

We can see that resources fluctuate around relatively stable equilibria. Resource consumption is low relative to the total concentrations, especially in the low-resource sites, so mostly resources are determined by transport.

## resource plot
plot(model_run, variable = "resources")

Here we can see the species distributions at the end of the simulation. Note that flume objects store the network in a list (to allow for duplicate networks to be run under the same conditions), hence the [[1]] index to plot the network.

## species distributions at the end
plot(model_run$network[[1]], variable = "species")

## additional plots -- forthcoming

## occupancy in time by reach
## this should be a map of the network, colored so that
## each node shows the proportion of time it was occupied by each species

## diversity over time
## diversity in space
## EF over time
## EF over space
## EF vs transport

flee-group/flume documentation built on Jan. 29, 2024, 6:44 p.m.