`EcologicalNetworksDynamics.Internals`

— ModuleEcologicalNetworksDynamics

Provide tools to simulate biomass dynamics in trophic and multiplex networks. Trophic networks only include feeding interactions, while multiplex networks also include non-trophic interactions such as interference between predators, or plant facilitation. The basic workflow has been designed to be as simple as possible, while remaining flexible for the experienced or adventurous user who would like to refine the model and its parameters.

Example of a simple workflow:

```
trophic_backbone = FoodWeb([1 => [2, 3]]) # species 1 eats plants 2 and 3
multi_net = MultiplexNetwork(trophic_backbone; L_facilitation = 1) # add 1 facilitation link
p = ModelParameters(multi_net) # generate model parameters
sol = simulate(p, rand(3)) # run simulation with random initial conditions
```

For more information, either go through the online documentation at (https://doc-url) or if you are looking for the help of a specific function read its docstring by writing `?<function_name>`

in a Julia REPL.

`EcologicalNetworksDynamics.Internals.AllometricParams`

— Type`AllometricParams(aₚ, aₑ, aᵢ, bₚ, bₑ, bᵢ)`

Parameters used to compute allometric rates for different metabolic classes.

The rate R is expressed as follow: $R = aMᵇ$, where a and b can take different values depending on the metabolic class of the species. This struct aims at storing these values of a and b. Specifically:

- aₚ: a for producers
- aₑ: a for ectotherm vertebrates
- aᵢ: a for invertebrates
- bₚ: b for producers
- bₑ: b for ectotherm vertebrates
- bᵢ: b for invertebrates

Default parameters values taken from the literature for certain rates can be accessed by calling the corresponding function, i.e. for:

- growth rate (r) call
`DefaultGrowthParams`

- metabolic rate (x) call
`DefaultMetabolismParams`

- max consumption rate (y) call
`DefaultMaxConsumptionParams`

**Example**

```
julia> params = AllometricParams(1, 2, 3, 4, 5, 6)
AllometricParams(aₚ=1, aₑ=2, aᵢ=3, bₚ=4, bₑ=5, bᵢ=6)
julia> params.aₚ
1
julia> params.aₑ
2
```

`EcologicalNetworksDynamics.Internals.BioRates`

— Method```
BioRates(
network::EcologicalNetwork;
d = allometric_rate(foodweb, DefaultMortalityParams()),
r = allometric_rate(foodweb, DefaultGrowthParams()),
x = allometric_rate(foodweb, DefaultMetabolismParams()),
y = allometric_rate(foodweb, DefaultMaxConsumptionParams()),
)
```

Compute the biological rates (r, x, y and e) of each species in the system.

The rates are:

- d: the natural mortality rate
- r: the growth rate
- x: the metabolic rate or metabolic demand
- y: the maximum consumption rate
- e: the assimilation efficiency If no value are provided for the rates, they take default values assuming an allometric scaling. Custom values can be provided for one or several rates by giving a vector of length 1 or S (species richness). Moreover, if one want to use allometric scaling ($R = aMᵇ$) but do not want to use default values for a and b, one can simply call
`allometric_rate`

with custom`AllometricParams`

.

**Examples**

```
julia> foodweb = FoodWeb([0 1; 0 0]); # sp. 1 "invertebrate", sp. 2 "producer"
julia> BioRates(foodweb) # default
BioRates:
d: [0.0, 0.0]
r: [0.0, 1.0]
x: [0.314, 0.0]
y: [8.0, 0.0]
e: (2, 2) sparse matrix
julia> BioRates(foodweb; r = [1.0, 1.0]) # specify custom vector for growth rate
BioRates:
d: [0.0, 0.0]
r: [1.0, 1.0]
x: [0.314, 0.0]
y: [8.0, 0.0]
e: (2, 2) sparse matrix
julia> BioRates(foodweb; x = 2.0) # if single value, fill the rate vector with it
BioRates:
d: [0.0, 0.0]
r: [0.0, 1.0]
x: [2.0, 2.0]
y: [8.0, 0.0]
e: (2, 2) sparse matrix
julia> custom_params = AllometricParams(3, 0, 0, 0, 0, 0); # use custom allometric params...
julia> BioRates(foodweb; y = allometric_rate(foodweb, custom_params)) # ...with allometric_rate
BioRates:
d: [0.0, 0.0]
r: [0.0, 1.0]
x: [0.314, 0.0]
y: [0.0, 3.0]
e: (2, 2) sparse matrix
```

`EcologicalNetworksDynamics.Internals.BioenergeticResponse`

— Method`BioenergeticResponse(B, i, j)`

Compute the bionergetic functional response for predator `i`

eating prey `j`

, given the species biomass `B`

. The bionergetic functional response is written:

\[F_{ij} = \frac{\omega_{ij} B_j^h}{B_0^h + c_i B_i B_0^h + \sum_k \omega_{ik} B_k^h}\]

With:

- $\omega$ the preferency, by default we assume that predators split their time equally among their preys, i.e. ∀j $\omega_{ij} = \omega_{i} = \frac{1}{n_{i,preys}}$ where $n_{i,preys}$ is the number of preys of predator i.
- $h$ the hill exponent, if $h = 1$ the functional response is of type II, and of type III if $h = 2$
- $c_i$ the intensity of predator intraspecific inteference
- $B_0$ the half-saturation density.

**Examples**

```
julia> foodweb = FoodWeb([0 0; 1 0]);
julia> F = BioenergeticResponse(foodweb)
BioenergeticResponse:
B0: [0.5, 0.5]
c: [0.0, 0.0]
h: 2.0
ω: (2, 2) sparse matrix
julia> F([1, 1], 1, 2) # no interaction, 1 does not eat 2
0.0
julia> F([1, 1], 2, 1) # interaction, 2 eats 1
0.8
julia> F([1.5, 1], 2, 1) # increases with resource biomass
0.9
```

See also `ClassicResponse`

, `LinearResponse`

and `FunctionalResponse`

.

`EcologicalNetworksDynamics.Internals.ClassicResponse`

— Method`ClassicResponse(B, i, j, mᵢ)`

Compute the classic functional response for predator `i`

eating prey `j`

, given the species biomass `B`

and the predator body mass `mᵢ`

. The classic functional response is written:

\[F_{ij} = \frac{1}{m_i} \cdot \frac{\omega_{ij} a_{r,ij} B_j^h} {1 + c_i B_i + h_t \sum_k \omega_{ik} a_{r,ik} B_k^h}\]

With:

- $\omega$ the preferency, by default we assume that predators split their time equally among their preys, i.e. ∀j $\omega_{ij} = \omega_{i} = \frac{1}{n_{i,preys}}$ where $n_{i,preys}$ is the number of preys of predator i.
- $h$ the hill exponent, if $h = 1$ the functional response is of type II, and of type III if $h = 2$
- $c_i$ the intensity of predator intraspecific inteference
- $a_{r,ij}$ the attack rate of predator i on prey j
- $h_t$ the handling time of predators
- $m_i$ the body mass of predator i

**Examples**

```
julia> foodweb = FoodWeb([0 0; 1 0]);
julia> F = ClassicResponse(foodweb; hₜ = 1.0, aᵣ = 0.5)
ClassicResponse:
c: [0.0, 0.0]
h: 2.0
ω: (2, 2) sparse matrix
hₜ: (2, 2) sparse matrix
aᵣ: (2, 2) sparse matrix
julia> F([1, 1], 1, 2, 1) # no interaction, 1 does not eat 2
0.0
julia> round(F([1, 1], 2, 1, 1); digits = 2) # interaction, 2 eats 1
0.33
julia> round(F([1.5, 1], 2, 1, 1); digits = 2) # increases with resource biomass
0.53
```

See also `BioenergeticResponse`

, `LinearResponse`

and `FunctionalResponse`

.

`EcologicalNetworksDynamics.Internals.Environment`

— Method`Environment(; T=293.15)`

Create environmental parameters of the system.

The environmental parameters are:

- T the temperature (in Kelvin) By default, the carrying capacities of producers are assumed to be 1 while capacities of consumers are assumed to be
`nothing`

as consumers do not have a growth term.

**Examples**

```
julia> e = Environment() # Default behaviour.
e.T == 293.15 # Kelvin
true
julia> e = Environment(; T = 300.0)
e.T == 300.0
true
```

See also `ModelParameters`

.

`EcologicalNetworksDynamics.Internals.ExponentialBA`

— Method`ExponentialBA()`

Generates the default parameters used to calculate temperature dependent rates using the exponential Boltzmann-Arrhenius method

`EcologicalNetworksDynamics.Internals.FoodWeb`

— Method```
FoodWeb(
A::AbstractMatrix;
species::Vector{<:Label} = default_speciesid(A),
Z::Real = 1,
M::Vector{<:Real} = compute_mass(A, Z),
metabolic_class::Vector{<:Label} = default_metabolic_class(A),
method::String = "unspecified",
)
```

**Generate a FoodWeb from an adjacency matrix**

This is the most direct way to generate a `FoodWeb`

. You only need to provide an adjacency matrix filled with 0s and 1s, that respectively indicates the absence or presence of an interaction between the corresponding species pair. For instance `A = [0 0; 1 0]`

corresponds to a system of 2 species in which species 2 eats species 1.

```
julia> FoodWeb([0 0; 1 0])
FoodWeb of 2 species:
A: sparse matrix with 1 links
M: [1.0, 1.0]
metabolic_class: 1 producers, 1 invertebrates, 0 vertebrates
method: unspecified
species: [s1, s2]
```

You can also provide optional arguments e.g. to change the species names.

```
julia> foodweb = FoodWeb([0 0; 1 0]; species = ["plant", "herbivore"]);
julia> foodweb.species == ["plant", "herbivore"]
true
```

**Generate a FoodWeb from an adjacency list**

An adjacency list is an iterable of `Pair`

s (e.g. vector of `Pair`

s) or a dictionary. If the adjacency list is an iterable of `Pair`

s, the first element of each pair is a predator and the second element of each pair are the preys eaten by the corresponding predator. If the adjacency list is a dictionary, keys are predators and values the corresponding preys.

Species can be identified either with `Integer`

s corresponding to species indexes or with labels (`String`

s or `Symbol`

s) corresponding to the species names. In the latter case, species will be ordered lexically. Moreover, if you use labels the species names will be directly passed to the `FoodWeb.species`

field.

```
julia> al_names = ["snake" => ("turtle", "mouse")]; # can also be `Symbol`s
julia> al_index = [2 => [1, 3]]; # ~ if sorting species lexically
julia> fw_from_names = FoodWeb(al_names);
julia> fw_from_index = FoodWeb(al_index);
julia> fw_from_names.A == fw_from_index.A == [0 0 0; 1 0 1; 0 0 0]
true
julia> fw_from_names.species == ["mouse", "snake", "turtle"] # ordered lexically
true
```

**Generate a FoodWeb from a structural model**

For larger communities it can be convenient to rely on a structural model. The structural models implemented are: `nichemodel(S; C)`

, `cascademodel(S; C)`

, nestedhierarchymodel(S; C) and `mpnmodel(S; C, p_forbidden)`

.

To generate a food web from one of these models you have to follow this syntax: `FoodWeb(model_name, model_arguments, optional_FoodWeb_arguments)`

.

For instance, to generate a `FoodWeb`

of 10 species (`S`

) with 15 links (`L`

) and predator-prey body mass ratio (`Z`

) of 50 with the niche model, you can do:

```
julia> foodweb = FoodWeb(nichemodel, 15; L = 15, Z = 50);
julia> richness(foodweb) == 15
true
julia> foodweb.method == "nichemodel"
true
```

Moreover, while generating the `FoodWeb`

we check that it does not have cycles or disconnected species. Moreover, while generating the `FoodWeb`

, by default we check that it does not disconnected species, but we do not check that it does not contain loops. These behaviours can respectively be changed by setting the keyword arguments `check_cycle`

and `check_disconnected`

to `true`

or `false`

.

By default, we set a tolerance corresponding to 10% of the number of links given in argument (`L`

), or 10% of the connectance if you give a connectance (`C`

). For instance, if you create a `FoodWeb`

with 20 links, by default we ensure that `FoodWeb`

in output will contain between 18 and 22 links.

```
julia> foodweb = FoodWeb(nichemodel, 15; L = 15, Z = 50); # by default tol_L = round(0.1*L)
julia> 13 <= n_links(foodweb) <= 18
true
```

However, the default tolerance can be changed with the `tol_L`

argument if you give a number of links (`L`

) or with `tol_C`

if you give a connectance (`C`

).

```
julia> foodweb = FoodWeb(nichemodel, 15; L = 15, Z = 50, tol_L = 0);
julia> n_links(foodweb) == 15
true
```

`FoodWeb`

s are generated randomly and can be rejected if they do not satisfy one of the properties defined above (number of links not in the tolerance range, cycles, etc.). Thus we repeat the generation of the `FoodWeb`

until all criteria are satisfied or until the number of maximum iteraction `iter_max`

is reached. If the maximum number of iterations is reached an error is thrown. That number can be controlled by the `iter_max`

keyword argument, default is set to `1e5`

.

**Generate a FoodWeb from a UnipartiteNetwork**

Lastly, EcologicalNetworkDynamics.jl has been designed to interact nicely with EcologicalNetworks.jl. Thus you can also create a `FoodWeb`

from a `UnipartiteNetwork`

:

```
julia> uni_net = cascademodel(10, 0.1); # generate a UnipartiteNetwork
julia> foodweb = FoodWeb(uni_net; quiet = true);
julia> foodweb.A == uni_net.edges # same adjacency matrices
true
```

`FoodWeb`

struct

The function returns a `FoodWeb`

which is a collection of the following fields:

`A`

the adjacency matrix`species`

the vector of species identities`M`

the vector of species body mass`metabolic_class`

the vector of species metabolic classes`method`

describes which model (e.g. niche model) was used to generate`A`

if no model has been used`method="unspecified"`

See also `MultiplexNetwork`

.

`EcologicalNetworksDynamics.Internals.FunctionalResponse`

— Method`FunctionalResponse(B)`

Compute functional response matrix given the species biomass `B`

. If `B`

is a scalar, all species are assumed to have the same biomass `B`

. Otherwise provide a vector s.t. `B[i]`

= biomass of species i.

**Examples**

```
julia> foodweb = FoodWeb([0 0; 1 0]);
julia> F = BioenergeticResponse(foodweb)
BioenergeticResponse:
B0: [0.5, 0.5]
c: [0.0, 0.0]
h: 2.0
ω: (2, 2) sparse matrix
julia> F([1, 1]) # provide a species biomass vector
2×2 SparseArrays.SparseMatrixCSC{Float64, Int64} with 1 stored entry:
⋅ ⋅
0.8 ⋅
julia> F(1) # or a scalar if homogeneous
2×2 SparseArrays.SparseMatrixCSC{Float64, Int64} with 1 stored entry:
⋅ ⋅
0.8 ⋅
julia> F([1.5, 1]) # response increases with resource biomass
2×2 SparseArrays.SparseMatrixCSC{Float64, Int64} with 1 stored entry:
⋅ ⋅
0.9 ⋅
```

See also `BioenergeticResponse`

, `LinearResponse`

and `ClassicResponse`

.

`EcologicalNetworksDynamics.Internals.GeneratedExpression`

— TypeWraps a Julia Expression generated by `generate_dbdt`

. Mostly useful for pretty-printing, but you can `eval`

uate it like a regular expression. The actual expression lies inside the `.expr`

attribute.

`EcologicalNetworksDynamics.Internals.Layer`

— Type`Layer(A::SparseMatrixCSC{Bool,Int64}, intensity::Union{Nothing,Float64})`

A `Layer`

of an interaction type contains two pieces of information concerning this interaction:

`A`

: where the interactions occur given by the adjacency matrix`intensity`

: the intensity of the interaction`f`

: the functional form of the non-trophic effect on the parameter associated with the layer

The intensity is only defined for non-trophic interactions and is set to `nothing`

for trophic interactions.

`EcologicalNetworksDynamics.Internals.LinearResponse`

— Method`LinearResponse(B, i, j)`

Compute the linear functional response for predator `i`

eating prey `j`

, given the species biomass `B`

. The linear functional response is written:

\[F_{ij} = \omega_{ij} \alpha_{i} B_j\]

With:

- $\omega$ the preferency, by default we assume that predators split their time equally among their preys, i.e. ∀j $\omega_{ij} = \omega_{i} = \frac{1}{n_{i,preys}}$ where $n_{i,preys}$ is the number of preys of predator i.
- $\alpha_{i}$ the consumption rate of predator i.

**Examples**

```
julia> foodweb = FoodWeb([0 0; 1 0]);
julia> F = LinearResponse(foodweb)
LinearResponse:
α: [⋅, 1.0]
ω: (2, 2) sparse matrix
julia> F([1, 1], 1, 2) # no interaction, 1 does not eat 2
0.0
julia> F([1, 1], 2, 1) # interaction, 2 eats 1
1.0
julia> F([1.5, 1], 2, 1) # increases linearly with resource biomass...
1.5
julia> F([1, 1.5], 2, 1) # but not with consumer biomass
1.0
```

See also `BioenergeticResponse`

, `ClassicResponse`

# Code generation version (raw) (↑ ↑ ↑ DUPLICATED FROM ABOVE ↑ ↑ ↑). and `FunctionalResponse`

.# (update together as long as the two coexist)

`EcologicalNetworksDynamics.Internals.LogisticGrowth`

— Type`LogisticGrowth`

constructor.

**Fields**

`a`

: producer competition matrix`K`

: vector of producer carrying capacities

See also `NutrientIntake`

.

`EcologicalNetworksDynamics.Internals.LogisticGrowth`

— Method```
LogisticGrowth(
network::EcologicalNetwork;
a = nothing,
K = 1,
)
```

Create the parameters for the logistic producer growth model. In the end, the `LogisticGrowth`

struct created stores a vector of carrying capacities `K`

and a competition matrix between the producers `a`

.

The carrying capacities can be specified via the `K`

arguments, by default they are all set to 1 for the producers and to `nothing`

for the consumers.

By default the competition matrix `a`

is assumed to have diagonal elements equal to 1 and 0 for all off-diagonal elements. This default can be changed via the `a`

argument. You can either directly pass the interaction matrix, a single `Number`

that will be interpreted as the value of all diagonal elements, a tuple of `Number`

the first will be interpreted as the diagonal value and the second as the off-diagonal values, or a named tuple of `Number`

. For the latter, the following aliases can be used to refer to the diagonal elements 'diagonal', 'diag', 'd', and for the off-diagonal elements 'offdiagonal', 'offdiag', 'o', 'rest', 'nondiagonal', 'nond'.

**Examples**

```
julia> foodweb = FoodWeb([0 0 0; 0 0 0; 1 1 0]) # 1 & 2 producers and 3 consumer.
g = LogisticGrowth(foodweb)
g.a == [
1 0 0
0 1 0
0 0 0
]
true
julia> g.K == [1, 1, nothing]
true
julia> K_cst = 5.0
g = LogisticGrowth(foodweb; K = K_cst) # Change the default K.
g.K == [K_cst, K_cst, nothing]
true
julia> K_vec = [1, 2, nothing]
g = LogisticGrowth(foodweb; K = K_vec) # Can also provide a vector.
g.K == K_vec
true
julia> g = LogisticGrowth(foodweb; a = (diag = 0.5, offdiag = 1.0))
g.a == [
0.5 1.0 0.0
1.0 0.5 0.0
0.0 0.0 0.0
]
true
julia> my_competition_matrix = [
0.5 1.0 0.0
1.0 0.5 0.0
0.0 0.0 0.0
]
g = LogisticGrowth(foodweb; a = my_competition_matrix)
g.a == my_competition_matrix
true
```

See also `ModelParameters`

and `NutrientIntake`

.

`EcologicalNetworksDynamics.Internals.ModelParameters`

— Method```
ModelParameters(
network::EcologicalNetwork;
biorates::BioRates=BioRates(foodweb),
environment::Environment=Environment(foodweb),
functional_response::FunctionalResponse=BioenergeticResponse(foodweb),
producer_growth::ProducerGrowth=ProducerGrowth(foodweb),
temperature_response::TemperatureResponse
)
```

Generate the parameters of the species community.

Default values are taken from Brose et al., 2006. The parameters are compartmented in different groups:

`FoodWeb`

: foodweb information (e.g. adjacency matrix)`BioRates`

: biological species rates (e.g. growth rates)`Environment`

: environmental variables (e.g. carrying capacities)`FunctionalResponse`

(F): functional response form (e.g. classic or bioenergetic functional response)`ProducerGrowth`

: model for producer growth (e.g. logistic or nutrient intake)`TemperatureResponse`

: method used for temperature dependency

**Examples**

```
julia> foodweb = FoodWeb([0 1; 0 0])
p = ModelParameters(foodweb)
ModelParameters{BioenergeticResponse, LogisticGrowth}:
network: FoodWeb(S=2, L=1)
environment: Environment(T=293.15K)
biorates: BioRates(d, r, x, y, e)
functional_response: BioenergeticResponse
producer_growth: LogisticGrowth
temperature_response: NoTemperatureResponse
julia> p.network # Check that stored foodweb is the same as the one we provided.
FoodWeb of 2 species:
A: sparse matrix with 1 links
M: [1.0, 1.0]
metabolic_class: 1 producers, 1 invertebrates, 0 vertebrates
method: unspecified
species: [s1, s2]
julia> p.functional_response # Default is bioenergetic.
BioenergeticResponse:
B0: [0.5, 0.5]
c: [0.0, 0.0]
h: 2.0
ω: (2, 2) sparse matrix
julia> classic_response = ClassicResponse(foodweb)
p = ModelParameters(foodweb; functional_response = classic_response);
p.functional_response # Check that the functional response is now "classic".
ClassicResponse:
c: [0.0, 0.0]
h: 2.0
ω: (2, 2) sparse matrix
hₜ: (2, 2) sparse matrix
aᵣ: (2, 2) sparse matrix
```

`ModelParameters`

can also be generated from a `MultiplexNetwork`

.

```
julia> foodweb = FoodWeb([0 1; 0 0])
net = MultiplexNetwork(foodweb)
p = ModelParameters(net; functional_response = ClassicResponse(net))
ModelParameters{ClassicResponse, LogisticGrowth}:
network: MultiplexNetwork(S=2, Lt=1, Lc=0, Lf=0, Li=0, Lr=0)
environment: Environment(T=293.15K)
biorates: BioRates(d, r, x, y, e)
functional_response: ClassicResponse
producer_growth: LogisticGrowth
temperature_response: NoTemperatureResponse
```

`EcologicalNetworksDynamics.Internals.MultiplexNetwork`

— Method`MultiplexNetwork(foodweb::FoodWeb; args...)`

Build the `MultiplexNetwork`

from a `FoodWeb`

. A multiplex is composed of a `trophic_layer`

and several non-trophic `Layer`

indexed by interaction types: (`:competition`

, `:facilitation`

, `:interference`

and `:refuge`

). The `:trophic`

layer is given by the `FoodWeb`

, and other, non-trophic layers are empty by default.

There are various ways to specify non-trophic layers to be not-empty, using variably-named `args`

whose general form is either:

```
<parameter_name>_<interaction_name> = value
eg. intensity_refuge = 4.0
or equivalently: I_r = 4.0
```

to set up `parameter`

for layer `interaction`

, or:

```
<parameter_name> = (<interaction1>=value1, <interaction2>=value2, ...)
eg. intensity = (facilitation = 2.0, interference = 1.5)
or equivalently: I = (f = 2.0, i = 1.5)
```

to set up `parameter`

for all layers `interaction1`

, `interaction2`

, etc., or again:

```
<interaction_name> = (<parameter1>=value1, <parameter2>=value2, ...)
eg. competition = (connectance = 4, symmetry = false)
or equivalently: c = (L = 4, s = false)
```

to set up `parameter1`

and `parameter2`

for the layer `interaction.

Valid interactions and parameters names can be retrieved with the following two cheat-sheets:

```
julia> interaction_names()
OrderedCollections.OrderedDict{Symbol, Vector{Symbol}} with 5 entries:
:trophic => [:t, :trh]
:competition => [:c, :cpt]
:facilitation => [:f, :fac]
:interference => [:i, :itf]
:refuge => [:r, :ref]
julia> multiplex_network_parameters_names()
OrderedCollections.OrderedDict{Symbol, Vector{Symbol}} with 6 entries:
:adjacency_matrix => [:A, :matrix, :adj_matrix]
:intensity => [:I, :int]
:functional_form => [:F, :fn]
:connectance => [:C, :conn]
:number_of_links => [:L, :n_links]
:symmetry => [:s, :sym, :symmetric]
```

**Fill non-trophic layers**

Specify a number of desired links to generate a non-trophic layer.

```
julia> foodweb = FoodWeb([0 0 0; 0 0 0; 1 1 0]); # 2 producers and 1 consumer
julia> MultiplexNetwork(foodweb; L_facilitation = 2) # or 'n_links_f', or 'L_fac', etc.
MultiplexNetwork of 3 species:
trophic_layer: 2 links
competition_layer: 0 links
facilitation_layer: 2 links
interference_layer: 0 links
refuge_layer: 0 links
```

Alternately, specify desired connectance for the layer (a value of `1.0`

creates as many links as possible).

```
julia> MultiplexNetwork(foodweb; C_cpt = 1.0) # or 'C_competition', or 'connectance_c', etc.
MultiplexNetwork of 3 species:
trophic_layer: 2 links
competition_layer: 2 links
facilitation_layer: 0 links
interference_layer: 0 links
refuge_layer: 0 links
```

Alternately, provide an explicit adjacency matrix.

```
julia> foodweb = FoodWeb([0 0 0; 1 0 0; 1 0 0]); # 2 consumers feeding on producer 1
julia> MultiplexNetwork(foodweb; A_interference = [0 0 0; 0 0 1; 0 1 0]) # or 'matrix_i' etc.
MultiplexNetwork of 3 species:
trophic_layer: 2 links
competition_layer: 0 links
facilitation_layer: 0 links
interference_layer: 2 links
refuge_layer: 0 links
```

**Set non-trophic intensity**

Non-trophic layers intensities default to '1.0'. Modify them with corresponding arguments:

```
julia> foodweb = FoodWeb([0 0 0; 0 0 0; 1 1 0]); # 2 producers and 1 consumer
julia> multiplex_network = MultiplexNetwork(foodweb; facilitation = (L = 1, I = 2.0));
julia> multiplex_network.layers[:facilitation] # or [:f] or [:fac] etc.
Layer(A=AdjacencyMatrix(L=1), intensity=2.0)
```

**Change assumptions about the symmetry of non-trophic interactions**

An interaction is symmetric iff "$i$ interacts with $j$" implies that "$j$ interacts with $i$". In other words, an interaction is symmetric iff the adjacency matrix of that interaction is symmetric.

With default settings:

- competition is symmetric
- facilitation is not symmetric
- interference is symmetric
- refuge is not symmetric

Change the defaults with the appropriate arguments.

For example, competition is assumed to be symmetric, then the number of competition links has to be even. But you can change this default as follows:

```
julia> foodweb = FoodWeb([0 0 0; 0 0 0; 1 1 0]);
julia> MultiplexNetwork(foodweb; competition = (sym = false, L = 1))
MultiplexNetwork of 3 species:
trophic_layer: 2 links
competition_layer: 1 links
facilitation_layer: 0 links
interference_layer: 0 links
refuge_layer: 0 links
```

If you don't specify `sym=false`

an error will be thrown.

**Construct arguments dynamically.**

Sometimes you need to dynamically choose the layer or parameters in a group. Use `Dict`

s instead of named tuples to this end:

```
julia> MultiplexNetwork(foodweb; f = (L = 1, I = 2.0))
MultiplexNetwork of 3 species:
trophic_layer: 2 links
competition_layer: 0 links
facilitation_layer: 1 links
interference_layer: 0 links
refuge_layer: 0 links
julia> dict = Dict(:L => 1, :I => 2.0);
b = MultiplexNetwork(foodweb; f = dict)
MultiplexNetwork of 3 species:
trophic_layer: 2 links
competition_layer: 0 links
facilitation_layer: 1 links
interference_layer: 0 links
refuge_layer: 0 links
```

`EcologicalNetworksDynamics.Internals.NutrientIntake`

— Method```
NutrientIntake(
network::EcologicalNetwork;
n_nutrients = 2,
turnover = 0.25,
supply = fill(10.0, n_nutrients),
concentration = collect(LinRange(1, 0.5, n_nutrients)),
half_saturation = 1.0,
)
```

Create parameters for the nutrient intake model of producer growths. These parameters are:

- the nutrient
`turnover`

rate (vector of length`n_nutrients`

) - the
`supply`

concentration for each nutrient (vector of length`n_nutrients`

) - the relative
`concentration`

of each nutrients in each producers (matrix of size (n*producers, n*nutrients)) - the
`half_saturation`

densities for each nutrient-producer pair (matrix of size (n*producers, n*nutrients))

**Examples**

```
julia> foodweb = FoodWeb([4 => [1, 2, 3]]); # 4 feeds on 1, 2 and 3.
ni = NutrientIntake(foodweb) # knights who say..
length(ni) # 2 nutrients by default.
2
julia> ni.turnover == [0.25, 0.25]
true
julia> ni.supply == [10, 10]
true
julia> ni.concentration == [1 0.5; 1 0.5; 1 0.5]
true
julia> ni.half_saturation == [1 1; 1 1; 1 1]
true
julia> ni = NutrientIntake(foodweb; turnover = 1.0)
ni.turnover == [1.0, 1.0]
true
julia> ni = NutrientIntake(foodweb; n_nutrients = 5)
length(ni)
5
```

See also `ModelParameters`

and `ProducerGrowth`

.

`Base.convert`

— MethodConvert a `MultiplexNetwork`

to a `FoodWeb`

. The convertion consists in removing the non-trophic layers of the multiplex network.

`Base.filter`

— MethodFilter species of the network (`net`

) for which `f(species_index, net) = true`

.

`Base.length`

— Method`length(n::NutrientIntake)`

Number of nutrients in the `n`

utrient intake model.

`Base.map`

— MethodTransform species of the network (`net`

) by applying `f`

to each species.

`Base.show`

— MethodOne line ModelParameters display.

`Base.show`

— MethodOne line `BioRates`

display.

`Base.show`

— MethodOne line Environment display.

`Base.show`

— MethodOne line FoodWeb display.

`Base.show`

— MethodOne line display FunctionalResponse

`Base.show`

— MethodOne line `Layer`

display.

`Base.show`

— MethodOne line `MultiplexNetwork`

display.

`Base.show`

— MethodOne line display FunctionalResponse

`Base.show`

— Method`Base.show(io::IO, temperature_response::TemperatureResponse)`

One-line TemperatureResponse display.

`Base.show`

— MethodMultiline ModelParameters display.

`Base.show`

— MethodMultiline `BioRates`

display.

`Base.show`

— MethodMultiline BioenergeticResponse display.

`Base.show`

— MethodMultiline ClassicResponse display.

`Base.show`

— MethodMultiline Environment display.

`Base.show`

— MethodMultiline TemperatureResponse::ExponentialBA display.

`Base.show`

— MethodMultiline FoodWeb display.

`Base.show`

— MethodMultiline LinearResponse display.

`Base.show`

— MethodMultiline LogisticGrowth display.

`Base.show`

— MethodMultiline `MultiplexNetwork`

display.

`Base.show`

— MethodMultiline NutrientIntake display.

`EcologicalNetworksDynamics.Internals.A_competition_full`

— Method`A_competition_full(net)`

Adjacency matrix of all possible competition links in the network `net`

.

**Example**

```
julia> foodweb = FoodWeb([0 0; 0 0]; quiet = true); # 2 producers
julia> A_competition_full(foodweb)
2×2 SparseArrays.SparseMatrixCSC{Float64, Int64} with 2 stored entries:
⋅ 1.0
1.0 ⋅
```

See also `A_interference_full`

, `A_facilitation_full`

, `A_refuge_full`

.

`EcologicalNetworksDynamics.Internals.A_facilitation_full`

— Method`A_facilitation_full(net)`

Adjacency matrix of all possible facilitation links in the network `net`

.

**Example**

```
julia> foodweb = FoodWeb([0 0; 0 0]; quiet = true); # 2 producers
julia> A_facilitation_full(foodweb)
2×2 SparseArrays.SparseMatrixCSC{Float64, Int64} with 2 stored entries:
⋅ 1.0
1.0 ⋅
```

See also `A_competition_full`

, `A_interference_full`

, `A_refuge_full`

.

`EcologicalNetworksDynamics.Internals.A_interference_full`

— Method`A_interference_full(net)`

Adjacency matrix of all possible interference links in the network `net`

.

**Example**

```
julia> foodweb = FoodWeb([0 0 0; 1 0 0; 1 0 0]); # 2 consumers eating producer 1
julia> A_interference_full(foodweb)
3×3 SparseArrays.SparseMatrixCSC{Float64, Int64} with 2 stored entries:
⋅ ⋅ ⋅
⋅ ⋅ 1.0
⋅ 1.0 ⋅
```

See also `A_competition_full`

, `A_facilitation_full`

, `A_refuge_full`

.

`EcologicalNetworksDynamics.Internals.A_nti_full`

— MethodAdjacency matrix of potential links given by the `potential_links`

function in `net`

.

`EcologicalNetworksDynamics.Internals.A_refuge_full`

— Method`A_refuge_full(net)`

Adjacency matrix of all possible refuge links in the network `net`

.

**Example**

```
julia> foodweb = FoodWeb([0 0 0; 0 0 0; 1 0 0]; quiet = true); # consumer 3 eats producer 1
julia> A_refuge_full(foodweb)
3×3 SparseArrays.SparseMatrixCSC{Float64, Int64} with 1 stored entry:
⋅ ⋅ ⋅
1.0 ⋅ ⋅
⋅ ⋅ ⋅
```

See also `A_competition_full`

, `A_interference_full`

, `A_facilitation_full`

.

`EcologicalNetworksDynamics.Internals.DefaultGrowthParams`

— Method`DefaultGrowthParams()`

Default allometric parameters (a, b) values for growth rate (r).

See also `AllometricParams`

`EcologicalNetworksDynamics.Internals.DefaultMaxConsumptionParams`

— Method`DefaultMaxConsumptionParams()`

Default allometric parameters (a, b) values for max consumption rate (y).

See also `AllometricParams`

`EcologicalNetworksDynamics.Internals.DefaultMetabolismParams`

— Method`DefaultMetabolismParams()`

Default allometric parameters (a, b) values for metabolic rate (x).

See also `AllometricParams`

`EcologicalNetworksDynamics.Internals.DefaultMortalityParams`

— Method`DefaultMortalityParams()`

Default allometric parameters (a, b) values for mortality rate (d).

See also `AllometricParams`

`EcologicalNetworksDynamics.Internals.ExtinctionCallback`

— Method`ExtinctionCallback(extinction_threshold::AbstractVector, verbose::Bool)`

Generate a DiffEqCallbacks.DiscreteCallback to extinguish species below the `extinction_threshold`

. The `extinction_threshold`

can be either: a `Number`

(same threshold for every species) or an `AbstractVector`

of length species richness (one threshold per species). If `verbose = true`

a message is printed when a species goes extinct, otherwise no message are printed.

`EcologicalNetworksDynamics.Internals.allometric_rate`

— MethodCompute rate vector (one value per species) with allometric scaling.

`EcologicalNetworksDynamics.Internals.allometricparams_to_vec`

— MethodCreate species parameter vectors for a, b of length S (species richness) given the allometric parameters for the different metabolic classes (aₚ,aᵢ,...).

`EcologicalNetworksDynamics.Internals.allometricscale`

— MethodAllometric scaling: parameter expressed as a power law of body-mass (M).

`EcologicalNetworksDynamics.Internals.asymmetrize`

— Method`asymmetrize(V)`

Remove duplicate tuples from a symmetric vector of tuples. A vector `V`

of tuples is said to be symmetric ⟺ ((i,j) ∈ `V`

⟺ (j,i) ∈ `V`

). The tuple that has its 1st element less than its 2nd element is kept i.e. if i < j (i,j) is kept, and (j,i) otherwise.

See also `symmetrize`

.

`EcologicalNetworksDynamics.Internals.attack_rate`

— Method`attack_rate(network::EcologicalNetwork)`

Compute the attack rate for all predator-prey couples of the system. The output `aᵣ`

is square matrix with length equal to the species richness with `aᵣ[i,j]`

corresponding to the attack rate of predator $i$ on prey $j$. The attack rate of a predator-prey couple is given by their body masses, formally:

- $a_{r,ij} = 50 m_i^{0.45} m_j^{0.15}$ if both species are mobiles;
- $a_{r,ij} = 50 m_j^{0.15}$ if i is sessile and j mobile;
- $a_{r,ij} = 50 m_i^{0.45}$ if j is sessile and i mobile.

This formula is taken from Miele et al. 2019 (PLOS Computational).

`EcologicalNetworksDynamics.Internals.biomass`

— Method`biomass(solution; kwargs...)`

Returns a named tuple of total and species biomass, averaged over the `last`

timesteps.

**Arguments**

`kwargs...`

arguments are forwarded to `extract_last_timesteps`

. See `extract_last_timesteps`

for the argument details.

Can also handle a species x time biomass matrix, e.g. `biomass(mat::AbstractMatrix;)`

or a vector, e.g. `biomass(vec::AbstractVector;)`

.

**Examples**

```
julia> foodweb = FoodWeb([0 0; 1 0]);
params = ModelParameters(foodweb);
B0 = [0.5, 0.5];
sol = simulate(params, B0);
bm = biomass(sol; last = 2);
biomass(sol; last = 2).species ≈ [0.1890006203352691, 0.21964742227673806]
true
julia> biomass(sol; last = 2, idxs = [1]); # Get biomass for one species only
[biomass(sol; last = 2, idxs = [1]).total] ≈
biomass(sol; last = 2, idxs = [1]).species ≈
[0.1890006203352691]
true
julia> biomass([2 1; 4 2])
(total = 4.5, species = [1.5, 3.0])
```

`EcologicalNetworksDynamics.Internals.boltzmann`

— Method`boltzmann(Ea, T, T0)`

Calculates the temperature dependence term of the Boltzmann-Arrhenius equation, normalised to 20°C.

`EcologicalNetworksDynamics.Internals.cascade_model`

— Method`cascade_model(S::Int, C::AbstractFloat)`

Generate an adjancency matrix using the cascade model from a number of species `S`

and a connectance `C`

.

**Examples**

`cascade_model(10, 0.2)`

See Cohen et al. (1985) for details.

`EcologicalNetworksDynamics.Internals.cascade_model`

— Method`cascade_model(S::Int, L::Int)`

Generate an adjancency matrix using the cascade model from a number of species `S`

and a number of links `L`

.

**Examples**

`cascade_model(10, 3)`

See Cohen et al. (1985) for details.

`EcologicalNetworksDynamics.Internals.clean_labels`

— MethodCheck that labels have the correct format and convert them to `String`

s if needed.

`EcologicalNetworksDynamics.Internals.clean_metabolic_class`

— MethodCheck that provided metabolic classes are valid.

`EcologicalNetworksDynamics.Internals.coefficient_of_variation`

— Method```
coefficient_of_variation(
solution::Solution;
threshold = 0,
last = "10%",
corrected = true,
kwargs...,
)
```

Computes the Coefficient of Variation (CV) of community biomass and its partition in average species CV and synchrony, following Thibault & Connolly (2013).

The function excludes dead species, i.e. species that have an average biomass below `threshold`

over the `last`

timesteps (see `living_species`

). It avoids division by 0 in the CV computation. `last`

= "10%" by default. We set `corrected = true`

by default (see `Statistics.std`

) which computes an unbiaised estimator of the variance.

See `richness`

for the argument details.

**Reference**

Thibaut, L. M., & Connolly, S. R. (2013). Understanding diversity–stability relationships : Towards a unified model of portfolio effects. Ecology Letters, 16(2), 140‑150. https://doi.org/10.1111/ele.12019

**Examples**

```
julia> foodweb = FoodWeb([0 1 1; 0 0 0; 0 0 0]); # Two producers and one consumer
params = ModelParameters(foodweb);
B0 = [0.5, 0.5, 0.5];
sol = simulate(params, B0; verbose = true);
key = (:community, :species_mean, :synchrony);
s = coefficient_of_variation(sol; last = 10)[key];
keys(s), round.(values(s); digits = 2)
((:community, :species_mean, :synchrony), (0.02, 0.06, 0.1))
julia> B0 = [0, 0.5, 0.5]; # Two producers
sol = simulate(params, B0; verbose = true);
s = coefficient_of_variation(sol; last = 10)[key];
keys(s), round.(values(s); digits = 2)
((:community, :species_mean, :synchrony), (0.15, 0.15, 1.0))
```

`EcologicalNetworksDynamics.Internals.community_cv`

— Method```
community_cv(
solution::Solution;
threshold = 0,
last = "10%",
corrected = true,
kwargs...,
)
```

Compute the temporal Coefficient of Variation of community biomass.

See `coefficient_of_variation`

for the argument details.

`EcologicalNetworksDynamics.Internals.connectance`

— MethodConnectance of network: number of links / (number of species)^2

`EcologicalNetworksDynamics.Internals.consumption`

— MethodCompute consumption terms of ODEs.

`EcologicalNetworksDynamics.Internals.cpad`

— FunctionPad text to the center of output. This one seems missing from Julia now :( https://github.com/JuliaLang/julia/pull/23187

```
julia> cpad(9, 3)
" 9 "
julia> cpad("test", 1)
"test"
julia> cpad("test", 6, '-'; left = false)
"-test-"
julia> cpad("test", 6, '-'; left = true)
"-test-"
julia> cpad("test", 7, '-'; left = false)
"-test--"
julia> cpad("test", 7, '-'; left = true)
"--test-"
```

`EcologicalNetworksDynamics.Internals.draw_asymmetric_links`

— Method`draw_asymmetric_links(potential_links, C::AbstractFloat)`

Draw randomly from the list of `potential_links`

such that the link connectance is `C`

. Links are drawn asymmetrically, i.e. $i$ interacts with $j$ ⇏ $j$ interacts with $i$.

`EcologicalNetworksDynamics.Internals.draw_asymmetric_links`

— Method`draw_asymmetric_links(potential_links, L::Integer)`

Draw randomly `L`

links from the list of `potential_links`

. Links are drawn asymmetrically, i.e. $i$ interacts with $j$ ⇏ $j$ interacts with $i$.

`EcologicalNetworksDynamics.Internals.draw_symmetric_links`

— Method`draw_symmetric_links(potential_links, C::AbstractFloat)`

Draw randomly from the list of `potential_links`

such that the link connectance is `C`

. Links are drawn symmetrically, i.e. $i$ interacts with $j$ ⇒ $j$ interacts with $i$.

`EcologicalNetworksDynamics.Internals.draw_symmetric_links`

— Method`draw_symmetric_links(potential_links, L::Integer)`

Draw randomly `L`

links from the list of `potential_links`

. Links are drawn symmetrically, i.e. $i$ interacts with $j$ ⇒ $j$ interacts with $i$.

`EcologicalNetworksDynamics.Internals.dudt!`

— Method`dudt!(du, u, p, _)`

Compute the species and nutrient (when relevant) abundance derivatives `du`

, given the abundances `u`

and the model `p`

. The last silent argument is the time at which is evaluated the derivatives and is a requirement of DifferentialEquations.

`EcologicalNetworksDynamics.Internals.effect_competition`

— MethodEffect of competition on the net growth rate.

`EcologicalNetworksDynamics.Internals.effect_facilitation`

— MethodEffect of facilitation on intrinsic growth rate.

`EcologicalNetworksDynamics.Internals.effect_refuge`

— MethodEffect of refuge on attack rate.

`EcologicalNetworksDynamics.Internals.efficiency`

— Method`efficiency(network; e_herbivore=0.45, e_carnivore=0.85)`

Create the assimilation efficiency matrix (`Efficiency`

). `Efficiency[i,j]`

is the assimation efficiency of predator i eating prey j. A perfect efficiency corresponds to an efficiency of 1. The efficiency depends on the metabolic class of the prey:

- if prey is producter, efficiency is
`e_herbivore`

- otherwise efficiency is
`e_carnivore`

Default values are taken from Miele et al. 2019 (PLOS Comp.).

`EcologicalNetworksDynamics.Internals.evenness`

— Method`evenness(solution; threshold = 0, kwargs...)`

Computes the average Pielou evenness, over the `last`

timesteps.

`kwargs...`

arguments are forwarded to `extract_last_timesteps`

. See `extract_last_timesteps`

for the argument details.

Can also handle a vector, e.g. `evenness(n::AbstractVector; threshold = 0)`

**Reference**

https://en.wikipedia.org/wiki/Species_evenness

`EcologicalNetworksDynamics.Internals.exp_ba_attack_rate`

— Method`exp_ba_attack_rate()`

Default temperature-dependent and allometric parameters (`a`

, `b`

, `c`

, `Eₐ`

) values for attack rate (aᵣ). (Rall et al., 2012, Binzer et al., 2016)

`EcologicalNetworksDynamics.Internals.exp_ba_carrying_capacity`

— Method`exp_ba_carrying_capacity()`

Default temperature-dependent and allometric parameters (`a`

, `b`

, `c`

, `Eₐ`

) values for carrying capacity. (Meehan, 2006, Binzer et al., 2016)

`EcologicalNetworksDynamics.Internals.exp_ba_growth`

— Method`exp_ba_growth()`

Default temperature-dependent and allometric parameters (`a`

, `b`

, `c`

, `Eₐ`

) values for growth rate (`r`

). (Savage et al., 2004, Binzer et al., 2016)

`EcologicalNetworksDynamics.Internals.exp_ba_handling_time`

— Method`exp_ba_handling_time()`

Default temperature-dependent and allometric parameters (`a`

, `b`

, `c`

, `Eₐ`

) values for handling time (`hₜ`

). (Rall et al., 2012, Binzer et al., 2016)

`EcologicalNetworksDynamics.Internals.exp_ba_matrix_rate`

— Method`exp_ba_matrix_rate(net, T, exponentialBAparams)`

Compute rate matrix (one value per species interaction) with temperature dependent scaling, given the parameters stored in an ExponentialBAParams struct `(aᵣ, hₜ)`

.

```
julia> fw = FoodWeb([0 0; 1 0]);
p = exp_ba_attack_rate();
exp_ba_matrix_rate(fw, 303.15, p)
2×2 SparseArrays.SparseMatrixCSC{Float64, Int64} with 1 stored entry:
⋅ ⋅
3.35932e-6 ⋅
```

`EcologicalNetworksDynamics.Internals.exp_ba_metabolism`

— Method`exp_ba_metabolism()`

Default temperature-dependent and allometric parameters (`a`

, `b`

, `c`

, `Eₐ`

) values for metabolic rate (`x`

). (Ehnes et al., 2011, Binzer et al., 2016)

`EcologicalNetworksDynamics.Internals.exp_ba_params_to_vec`

— Method`exp_ba_params_to_vec(foodweb, ExponentialBAparams)`

Create species parameter vectors for `a`

, `b`

, `c`

of length `S`

(species richness) given the parameters for the different metabolic classes (`aₚ`

, `a`

, ...).

`EcologicalNetworksDynamics.Internals.exp_ba_vector_rate`

— Method`exp_ba_vector_rate(net, T ExponentialBAParams)`

Compute rate vector (one value per species) with temperature dependent scaling, given the parameters stored in an ExponentialBAParams struct `(x, r, K)`

.

```
julia> fw = FoodWeb([0 0; 1 0])
p = exp_ba_growth()
exp_ba_vector_rate(fw, 303.15, p) ≈ [4.6414071636920824e-7, 0.0]
true
```

`EcologicalNetworksDynamics.Internals.extract_last_timesteps`

— Method`extract_last_timesteps(solution; idxs = nothing, quiet = false, kwargs...)`

Returns the biomass matrix of species x time over the `last`

timesteps.

**Arguments**

`last`

: the number of last timesteps to consider. A percentage can also be also be provided as a`String`

ending by`%`

. Defaulted to 1.`idxs`

: vector of species indexes or names. Set to`nothing`

by default.`quiet`

: ignores warning issue while extracting timesteps before last species extinction

See `richness`

for the other arguments. If `idxs`

is an integer, it returns a vector of the species biomass instead of a matrix.

**Examples**

```
julia> fw = FoodWeb([0 0; 1 0]);
p = ModelParameters(fw);
m = simulate(p, [0.5, 0.5]);
sim = extract_last_timesteps(m; last = 1, idxs = [2, 1]);
sim ==
extract_last_timesteps(m; last = 1, idxs = ["s2", "s1"]) ≈
[0.219659439; 0.188980349;;]
true
julia> sim1 = extract_last_timesteps(m; last = 1, idxs = [1, 2]);
sim2 = extract_last_timesteps(m; last = 1, idxs = ["s1", "s2"]);
sim1 ≈ sim2 ≈ [0.188980349; 0.219659439;;]
true
julia> sim = extract_last_timesteps(m; last = 1, idxs = [2]);
sim == extract_last_timesteps(m; last = 1, idxs = "s2")
true
```

`EcologicalNetworksDynamics.Internals.fill_sparsematrix`

— Method`fill_sparsematrix(scalar, template_matrix)`

Return a matrix filled with a constant (`scalar`

) for indexes where the value of the `template_matrix`

is non-zero.

**Examples**

```
julia> template_matrix = ones(2, 2);
julia> Internals.fill_sparsematrix(10, template_matrix)
2×2 SparseArrays.SparseMatrixCSC{Float64, Int64} with 4 stored entries:
10.0 10.0
10.0 10.0
julia> template_matrix[1, 1] = 0;
julia> Internals.fill_sparsematrix(10, template_matrix)
2×2 SparseArrays.SparseMatrixCSC{Float64, Int64} with 3 stored entries:
⋅ 10.0
10.0 10.0
```

`EcologicalNetworksDynamics.Internals.find_steady_state`

— Method`find_steady_state(params::ModelParameters, B0::AbstractVector; kwargs...)`

Find a steady state of a system parametrized by a parameter set `params`

and some initial conditions `B0`

. The function returns a tuple (`steady_state`

, `terminated`

). If no steady state has been found `terminated`

is set to `false`

and `steady_state`

to `nothing`

. Otherwise, `terminated`

is set to `true`

and `steady_state`

to the corresponding value.

If you are not only interested in the steady state biomass, but also in the trajectories see `simulate`

.

**Examples**

```
julia> foodweb = FoodWeb([0 0; 1 0]); # create the foodweb
julia> biorates = BioRates(foodweb; d = 0); # set natural death to 0
julia> params = ModelParameters(foodweb; biorates = biorates);
julia> B0 = [0.5, 0.5]; # set initial biomass
julia> solution = find_steady_state(params, B0);
julia> solution.terminated # => a steady state has been found
true
julia> round.(solution.steady_state, digits = 2) # steady state biomass
2-element Vector{Float64}:
0.19
0.22
```

`EcologicalNetworksDynamics.Internals.generate_dbdt`

— Method`generate_dbdt(parms::ModelParameters, type)`

Produce a specialized julia expression and associated data, supposed to improve efficiency of subsequent simulations. The returned expression is typically `eval`

uated then passed along with the data as a `diff_code_data`

argument to `simulate`

.

There are two possible code generation styles:

With

`type = :raw`

, the generated expression is a straightforward translation of the underlying differential equations, with no loops, no recursive calls, nor heap-allocations: only local variables and basic arithmetic is used. This makes simulation very efficient, but the length of the generated expression varies with the number of species interactions. When the length becomes high, it takes much longer for julia to compile it. If it takes forever, (typically over`SyntaxtTree.callcount(expression) > 20_000`

), wait until julia 1.9 (maybe) or use the alternate style instead.With

`type = :compact`

, the generated expression is a more sophisticated implementation of the underlying differential equations, involving carefully crafted minimal loops and exactly one fixed-size heap-allocated bunch of data, reused on every call during simulation. This makes the simulation slightly less efficient than the above but, as the expression size no longer depends on the number of species interactions, there is no limit to using it and speedup simulations.

`EcologicalNetworksDynamics.Internals.get_alive_species`

— Method`get_alive_species(solution; idxs = nothing, threshold = 0)`

Returns a tuple with species having a biomass above `threshold`

at the end of a simulation.

**Examples**

```
julia> foodweb = FoodWeb([0 0; 0 0]; quiet = true);
params = ModelParameters(foodweb);
sim = simulate(params, [0, 0.5]; tmax = 20);
get_alive_species(sim)
(species = ["s2"], idxs = [2])
julia> sim = simulate(params, [0.5, 0]; tmax = 20);
get_alive_species(sim)
(species = ["s1"], idxs = [1])
julia> sim = simulate(params, [0.5, 0.5]; tmax = 20);
get_alive_species(sim)
(species = ["s1", "s2"], idxs = [1, 2])
julia> sim = simulate(params, [0, 0]; tmax = 20);
get_alive_species(sim)
(species = String[], idxs = Int64[])
```

`EcologicalNetworksDynamics.Internals.get_extinct_species`

— Method`get_extinct_species(sol)`

Extract list of extinct species from the solution returned by `simulate()`

.

```
julia> foodweb = FoodWeb([0 0; 1 0]);
params = ModelParameters(foodweb);
sol = simulate(params, [0.5, 0.5]);
isempty(get_extinct_species(sol)) # no species extinct
true
```

See also `simulate`

, `get_parameters`

.

`EcologicalNetworksDynamics.Internals.get_parameters`

— Method`get_parameters(sol)`

Extract the `ModelParameters`

input from the solution returned by `simulate()`

.

```
julia> foodweb = FoodWeb([0 0; 1 0]);
params = ModelParameters(foodweb);
sol = simulate(params, [0.5, 0.5]);
isa(get_parameters(sol), ModelParameters)
true
```

See also `simulate`

, `get_extinct_species`

.

`EcologicalNetworksDynamics.Internals.get_trophic_adjacency`

— MethodReturn the adjacency matrix of the trophic interactions.

`EcologicalNetworksDynamics.Internals.handling_time`

— Method`handling_time(network::EcologicalNetwork)`

Compute the handling time for all predator-prey couples of the system. The output `hₜ`

is a square matrix with length equal to the species richness, with `hₜ[i,j]`

corresponding to the handling time of predator $i$ on prey $j$. The handling time of a predator-prey couple is given by their body masses, formally: $h_{t,ij} = 0.3 m_i^{-0.48} m_j^{-0.66}$. This formula is taken from Miele et al. 2019 (PLOS Computational) and Rall et al. 2012 (Phil. Tran. R. Soc. B).

`EcologicalNetworksDynamics.Internals.homogeneous_preference`

— Method`homogeneous_preference(network::EcologicalNetwork)`

Create the preferency matrix (`ω`

) which describes how each predator splits its time between its different preys. `ω[i,j]`

is the fraction of time of predator i spent on prey j. By definition, ∀i $\sum_j \omega_{ij} = 1$. Here we assume an **homogeneous** preference, meaning that each predator splits its time equally between its preys, i.e. ∀j $\omega_{ij} = \omega_{i} = \frac{1}{n_{preys,i}}$ where $n_{preys,i}$ is the number of prey of predator i.

`EcologicalNetworksDynamics.Internals.invertebrates`

— MethodReturn indices of the invertebrates of the network (`net`

).

`EcologicalNetworksDynamics.Internals.is_boostable`

— Method`is_boostable(parms::ModelParameters, type)`

Return true if boost has been implemented for this parametrization with the given boosting style.

`EcologicalNetworksDynamics.Internals.is_model_net_valid`

— MethodCheck that `net`

does not contain cycles and does not have disconnected nodes.

`EcologicalNetworksDynamics.Internals.isinvertebrate`

— MethodIs species `i`

an invertebrate?

`EcologicalNetworksDynamics.Internals.ispredator`

— MethodIs species `i`

of the network (`net`

) a predator?

`EcologicalNetworksDynamics.Internals.isprey`

— MethodIs species `i`

a prey?

`EcologicalNetworksDynamics.Internals.isproducer`

— MethodIs species `i`

of the network (`net`

) a producer?

`EcologicalNetworksDynamics.Internals.isvertebrate`

— MethodIs species `i`

an ectotherm vertebrate?

`EcologicalNetworksDynamics.Internals.living_species`

— Method`living_species(solution::Solution; threshold = 0, idxs = nothing, kwargs...)`

Returns the vectors of alive species and their indices in the original network. Living species are the ones having, in average, a biomass above `threshold`

over the `last`

timesteps. `kwargs...`

are optional arguments passed to `extract_last_timesteps`

.

`kwargs...`

arguments are forwarded to `extract_last_timesteps`

. See `extract_last_timesteps`

for the argument details.

**Examples**

```
julia> foodweb = FoodWeb([0 1 1; 0 0 0; 0 0 0]);
params = ModelParameters(foodweb);
B0 = [0.5, 0.5, 0.5];
sol = simulate(params, B0; verbose = true);
julia> living_species(sol)
(species = ["s1", "s2", "s3"], idxs = [1, 2, 3])
julia> B0 = [0, 0.5, 0.5];
sol = simulate(params, B0; verbose = true);
living_species(sol; last = 1)
(species = ["s2", "s3"], idxs = [2, 3])
```

`EcologicalNetworksDynamics.Internals.mass_ratio`

— Method`mass_ratio(p::ModelParameters)`

Mean predator-prey body mass ratio given the model `p`

arameters.

`EcologicalNetworksDynamics.Internals.mass_ratio`

— Method`mass_ratio(network::EcologicalNetwork)`

Mean predator-prey body mass ratio given the `network`

.

`EcologicalNetworksDynamics.Internals.max_trophic_level`

— Method```
max_trophic_level(solution::Solution; threshold = 0, kwargs...)
mean_trophic_level(solution::Solution; threshold = 0, kwargs...)
weighted_mean_trophic_level(solution::Solution; threshold = 0, kwargs...)
```

Return the aggregated trophic level over the `last`

timesteps, either with `max`

or `mean`

aggregation, or by the mean trophic level weighted by species biomasses.

kwargs... arguments are forwarded to `extract_last_timesteps`

. See `extract_last_timesteps`

for the argument details.

These functions also handle biomass vectors associated with a network, as well as a vector of trophic levels (See examples).

**Examples**

```
julia> A = [0 0; 1 0];
max_trophic_level(A)
2.0
julia> mean_trophic_level(A)
1.5
julia> bm = [1, 1];
max_trophic_level(bm, A)
2.0
julia> mean_trophic_level(bm, A)
1.5
julia> weighted_mean_trophic_level(bm, A)
1.5
julia> bm = [0, 1];
max_trophic_level(bm, A)
1.0
julia> mean_trophic_level(bm, A)
1.0
julia> weighted_mean_trophic_level(bm, A)
1.0
julia> foodweb = FoodWeb(A; quiet = true);
params = ModelParameters(foodweb);
B0 = [0.5, 0.5];
sol = simulate(params, B0; verbose = false);
max_trophic_level(sol)
2.0
julia> mean_trophic_level(sol)
1.5
julia> w = weighted_mean_trophic_level(sol);
round(w; digits = 2)
1.54
```

`EcologicalNetworksDynamics.Internals.mean_trophic_level`

— Method```
max_trophic_level(solution::Solution; threshold = 0, kwargs...)
mean_trophic_level(solution::Solution; threshold = 0, kwargs...)
weighted_mean_trophic_level(solution::Solution; threshold = 0, kwargs...)
```

Return the aggregated trophic level over the `last`

timesteps, either with `max`

or `mean`

aggregation, or by the mean trophic level weighted by species biomasses.

kwargs... arguments are forwarded to `extract_last_timesteps`

. See `extract_last_timesteps`

for the argument details.

These functions also handle biomass vectors associated with a network, as well as a vector of trophic levels (See examples).

**Examples**

```
julia> A = [0 0; 1 0];
max_trophic_level(A)
2.0
julia> mean_trophic_level(A)
1.5
julia> bm = [1, 1];
max_trophic_level(bm, A)
2.0
julia> mean_trophic_level(bm, A)
1.5
julia> weighted_mean_trophic_level(bm, A)
1.5
julia> bm = [0, 1];
max_trophic_level(bm, A)
1.0
julia> mean_trophic_level(bm, A)
1.0
julia> weighted_mean_trophic_level(bm, A)
1.0
julia> foodweb = FoodWeb(A; quiet = true);
params = ModelParameters(foodweb);
B0 = [0.5, 0.5];
sol = simulate(params, B0; verbose = false);
max_trophic_level(sol)
2.0
julia> mean_trophic_level(sol)
1.5
julia> w = weighted_mean_trophic_level(sol);
round(w; digits = 2)
1.54
```

`EcologicalNetworksDynamics.Internals.min_max`

— Method`min_max(solution; kwargs...)`

Returns the vectors of minimum and maximum biomass of each species over the `last`

timesteps.

`kwargs...`

arguments are forwarded to `extract_last_timesteps`

. See `extract_last_timesteps`

for the argument details.

**Examples**

```
julia> foodweb = FoodWeb([0 0; 1 0]);
params = ModelParameters(foodweb);
sol = simulate(params, [0.5; 0.5]; tmax = 100);
ti = min_max(sol; last = 10);
round.(ti.min, digits = 3) # Min
2-element Vector{Float64}:
0.168
0.206
julia> round.(ti.max, digits = 3) # Max
2-element Vector{Float64}:
0.196
0.221
```

`EcologicalNetworksDynamics.Internals.model_foodweb`

— MethodGenerate a food web of `S`

species and number of links `L`

from a structural `model`

. Loop until the generated has a number of links in [L - ΔL; L + ΔL]. If the maximum number of iterations is reached an error is thrown instead.

`EcologicalNetworksDynamics.Internals.model_foodweb_from_C`

— MethodGenerate a food web of `S`

species and connectance `C`

from a structural `model`

. Loop until the generated has connectance in [C - ΔC; C + ΔC].

`EcologicalNetworksDynamics.Internals.n_links`

— Method`n_links(network)`

Compute the number of links of a network. If the argument is an adjacency matrix or a `Layer`

, then an integer corresponding to the number of links (i.e. 1s) is returned. If the argument is a `FoodWeb`

, then an integer corresponding to the number of trophic links is returned. If the argument is a `MultiplexNetwork`

, then a dictionnary is returned where `Dict[:interaction_type]`

contains the number of links of the selected interactions

**Examples**

```
julia> foodweb = FoodWeb([0 0 0; 1 0 0; 0 1 0]); # food chain of length 3 (2 links)
julia> n_links(foodweb) == n_links(foodweb.A) == 2
true
julia> multi_net = MultiplexNetwork(foodweb; L_facilitation = 1); # + 1 facilitation link
julia> n_links(multi_net)
EcologicalNetworksDynamics.Internals.InteractionDict{Int64} with 5 entries:
:trophic => 2
:facilitation => 1
:competition => 0
:refuge => 0
:interference => 0
```

`EcologicalNetworksDynamics.Internals.niche_model`

— Method`niche_model(S::Int, C::AbstractFloat)`

Generate an adjancency matrix using the niche model from a number of species `S`

and a connectance `C`

.

**Example**

`niche_model(10, 0.2)`

See Williams and Martinez (2000) for details.

`EcologicalNetworksDynamics.Internals.niche_model`

— Method`niche_model(S::Int, C::AbstractFloat)`

Generate an adjancency matrix using the niche model from a number of species `S`

and a number of links `L`

.

**Example**

`niche_model(10, 20) # 20 links for 10 species.`

```
See [Williams and Martinez (2000)](https://doi.org/10.1038/35004572) for details.
```

`EcologicalNetworksDynamics.Internals.nontrophic_adjacency_matrix`

— Method```
nontrophic_adjacency_matrix(
foodweb::FoodWeb,
find_potential_links::Function,
n;
symmetric=false)
```

Generate a non-trophic adjacency matrix given the function finding potential links (`find_potential_links`

) and `n`

which is interpreted as a number of links if `n:<Integer`

or as a connectance if `n:<AbstractFloat`

.

See also `potential_competition_links`

, `potential_facilitation_links`

, `potential_interference_links`

, `potential_refuge_links`

.

`EcologicalNetworksDynamics.Internals.number_of_resource`

— MethodNumber of resources species `i`

is feeding on.

`EcologicalNetworksDynamics.Internals.number_of_resource`

— MethodReturn a vector where element i is the number of resource(s) of species i.

`EcologicalNetworksDynamics.Internals.nutrient_dynamics`

— Method`nutrient_dynamics(model::ModelParameters, B, i_nutrients, n, G)`

Compute the dynamics of the nutrient `i_nutrient`

given its abundance `n`

, the species biomass `B`

and the vector of species growths `G`

and the model `p`

.

The nutrient dynamics is applicable only if `p`

is of type `NutrientIntake`

.

`EcologicalNetworksDynamics.Internals.nutrient_indices`

— MethodNutrient indices of the given `network`

.

`EcologicalNetworksDynamics.Internals.nutrient_richness`

— Method`nutrient_richness(p::ModelParameters)`

Number of nutrients in the model `p`

.

```
julia> foodweb = FoodWeb([1 => 2]) # 1 eats 2.
producer_growth = LogisticGrowth(foodweb) # No nutrients.
model = ModelParameters(foodweb; producer_growth)
nutrient_richness(model)
0
julia> producer_growth = NutrientIntake(foodweb; n_nutrients = 3) # Include nutrients.
model_with_nutrients = ModelParameters(foodweb; producer_growth)
nutrient_richness(model_with_nutrients)
3
```

See also `richness`

and `total_richness`

.

`EcologicalNetworksDynamics.Internals.parse_pair`

— MethodParse pairs within `FoodWeb()`

method working on adjacency list.

`EcologicalNetworksDynamics.Internals.potential_competition_links`

— Method`potential_competition_links(foodweb::FoodWeb)`

Find possible competition links. Competition links can only occur between two sessile species (producers).

`EcologicalNetworksDynamics.Internals.potential_facilitation_links`

— Method`potential_facilitation_links(foodweb::FoodWeb)`

Find possible facilitation links. Facilitation links can only occur from any species to a plant.

`EcologicalNetworksDynamics.Internals.potential_interference_links`

— Method`potential_interference_links(foodweb::FoodWeb)`

Find possible interference links. Interference links can only occur between two predators sharing at least one prey.

`EcologicalNetworksDynamics.Internals.potential_refuge_links`

— Method`potential_refuge_links(foodweb::FoodWeb)`

Find possible refuge links. Refuge links can only occur from a sessile species (producer) to a prey.

`EcologicalNetworksDynamics.Internals.predators`

— MethodReturn indices of the predators of the given `network`

.

`EcologicalNetworksDynamics.Internals.predators_of`

— Method`predators_of(i, network)`

Return indices of the predators of species `i`

for the given `network`

.

**Examples**

```
julia> foodweb = FoodWeb([0 0 0; 1 0 0; 1 1 0]);
julia> predators_of(1, foodweb)
2-element Vector{Int64}:
2
3
julia> predators_of(2, foodweb)
1-element Vector{Int64}:
3
julia> predators_of(3, foodweb)
Int64[]
```

`EcologicalNetworksDynamics.Internals.preys`

— MethodReturn indices of the preys of the network (`net`

).

`EcologicalNetworksDynamics.Internals.preys_of`

— Method`preys_of(i, network)`

Return indices of the preys of species `i`

for the given `network`

.

**Examples**

```
julia> foodweb = FoodWeb([0 0 0; 0 0 0; 1 1 0]);
julia> preys_of(3, foodweb)
2-element Vector{Int64}:
1
2
julia> preys_of(1, foodweb) # empty
Int64[]
```

See also `predators_of`

and `producers`

.

`EcologicalNetworksDynamics.Internals.process_idxs`

— Method`process_idxs(solution; idxs = nothing)`

Check and sanitize the species indices or names provided (`idxs`

). Used in `extract_last_timesteps`

and `living_species`

.

`EcologicalNetworksDynamics.Internals.producer_growth`

— Method`producer_growth(solution; kwargs...)`

Returns the average growth rates of producers over the `last`

timesteps as as well as the average (`mean`

) and the standard deviation (`std`

). It also returns by all the growth rates (`all`

) as a species x timestep matrix (as the solution matrix).

kwargs... arguments are forwarded to `extract_last_timesteps`

. See `extract_last_timesteps`

for the argument details.

**Examples**

```
julia> foodweb = FoodWeb([0 1 1; 0 0 0; 0 0 0]);
params = ModelParameters(foodweb);
B0 = [0.5, 0.5, 0.5];
sol = simulate(params, B0);
g = producer_growth(sol; last = 10);
g.species, round.(g.mean, digits = 2)
(["s2", "s3"], [0.15, 0.15])
```

`EcologicalNetworksDynamics.Internals.producers`

— MethodReturn indices of the producers of the given `network`

.

`EcologicalNetworksDynamics.Internals.remove_species`

— Method`remove_species(net::EcologicalNetwork, species_to_remove)`

Remove a list of species from a `FoodWeb`

, a `MultiplexNetwork`

or an adjacency matrix.

**Examples**

Adjacency matrix.

```
julia> A = [0 0 0; 1 0 0; 1 0 0]
remove_species(A, [1]) == [0 0; 0 0]
true
julia> remove_species(A, [2]) == [0 0; 1 0]
true
```

`FoodWeb`

.

```
julia> foodweb = FoodWeb([0 0 0; 1 0 0; 1 0 0]; M = [1, 10, 20])
new_foodweb = remove_species(foodweb, [2])
new_foodweb.A == [0 0; 1 0]
true
julia> new_foodweb.M == [1, 20]
true
```

`MultiplexNetwork`

.

```
julia> foodweb = FoodWeb([0 0 0; 1 0 0; 1 0 0]; M = [1, 10, 20])
net = MultiplexNetwork(foodweb; A_facilitation = [0 0 0; 1 0 0; 0 0 0])
new_net = remove_species(net, [2])
new_net.layers[:trophic].A == [0 0; 1 0]
true
julia> new_net.layers[:facilitation].A == [0 0; 0 0]
true
julia> new_net.M == [1, 20]
true
```

`EcologicalNetworksDynamics.Internals.replace`

— MethodConstruct a copy of the expression with the replacements given in `rep`

.

```
julia> import EcologicalNetworksDynamics.Internals: replace
julia> replace(:(a + (b + c / a)), Dict(:a => 5, :b => 8))
:(5 + (8 + c / 5))
```

`EcologicalNetworksDynamics.Internals.replace_vertebrates!`

— MethodDoes the user want to replace 'vertebrates' by 'ectotherm vertebrates'?

`EcologicalNetworksDynamics.Internals.richness`

— Method`richness(p::ModelParameters)`

Number of species in the network.

```
julia> foodweb = FoodWeb([1 => 2]) # 1 eats 2.
model = ModelParameters(foodweb)
richness(model)
2
```

See also `nutrient_richness`

and `total_richness`

.

`EcologicalNetworksDynamics.Internals.richness`

— Method`richness(n::AbstractVector; threshold = 0)`

When applied to a vector of biomass, returns the number of biomass above `threshold`

**Examples**

```
julia> richness([0, 1])
1
julia> richness([1, 1])
2
```

`EcologicalNetworksDynamics.Internals.richness`

— Method`richness(solution::Solution; threshold = 0, kwargs...)`

Returns the average number of species with a biomass larger than `threshold`

over the `last`

timesteps. `kwargs...`

are optional arguments passed to `extract_last_timesteps`

.

**Arguments:**

`solution`

: output of`simulate()`

or`solve()`

`threshold`

: biomass threshold below which a species is considered extinct. Set to 0 by debault and it is recommended to let as this. It is recommended to change the threshold using`ExtinctionCallback`

in`simulate`

.

**Examples**

```
julia> foodweb = FoodWeb([0 0; 1 0]);
params = ModelParameters(foodweb);
B0 = [0.5, 0.5];
sol = simulate(params, B0);
richness(sol; last = 10)
2.0
julia> sha = shannon_diversity(sol);
round(sha; digits = 3)
0.69
julia> simp = simpson(sol);
round(simp; digits = 3)
0.353
julia> even = evenness(sol);
round(even; digits = 3)
0.996
```

`EcologicalNetworksDynamics.Internals.set_temperature!`

— Method`set_temperature!(p::ModelParameters, T, F!::TemperatureResponse)`

Parametrize the given model given a temperature value and the associated temperature_response.

```
julia> fw = FoodWeb([0 0; 1 0]);
model = ModelParameters(fw; functional_response = ClassicResponse(fw));
model.temperature_response
NoTemperatureResponse
julia> model.biorates # Default Biorates.
BioRates:
d: [0.0, 0.0]
r: [1.0, 0.0]
x: [0.0, 0.314]
y: [0.0, 8.0]
e: (2, 2) sparse matrix
julia> set_temperature!(model, 25, ExponentialBA());
model.temperature_response
Parameters for ExponentialBA response:
r: ExponentialBAParams(1.5497531357028967e-7, 0.0, 0.0, -0.25, -0.25, -0.25, 0.0, 0.0, 0.0, -0.84)
x: ExponentialBAParams(0.0, 6.557967639824989e-8, 6.557967639824989e-8, -0.31, -0.31, -0.31, 0.0, 0.0, 0.0, -0.69)
aᵣ: ExponentialBAParams(0.0, 2.0452306245234897e-6, 2.0452306245234897e-6, 0.25, 0.25, 0.25, -0.8, -0.8, -0.8, -0.38)
hₜ: ExponentialBAParams(0.0, 15677.784668089162, 15677.784668089162, -0.45, -0.45, -0.45, 0.47, 0.47, 0.47, 0.26)
K: ExponentialBAParams(3.0, nothing, nothing, 0.28, 0.28, 0.28, 0.0, 0.0, 0.0, 0.71)
julia> model.biorates # Biorates given this temperature response.
BioRates:
d: [0.0, 0.0]
r: [1.9446555905910437e-162, 0.0]
x: [0.0, 3.7697909572935493e-135]
y: [0.0, 8.0]
e: (2, 2) sparse matrix
```

`EcologicalNetworksDynamics.Internals.shannon_diversity`

— Method`shannon_diversity(solution; threshold = 0, kwargs...)`

Computes the average Shannon entropy index, i.e. the first Hill number, over the `last`

timesteps.

`kwargs...`

arguments are forwarded to `extract_last_timesteps`

. See `extract_last_timesteps`

for the argument details.

Can also handle a vector, e.g. shannon_diversity(n::AbstractVector; threshold = 0)

**Reference**

https://en.wikipedia.org/wiki/Diversity*index#Shannon*index

`EcologicalNetworksDynamics.Internals.share_prey`

— MethodDo species `i`

and `j`

share at least one prey?

`EcologicalNetworksDynamics.Internals.simpson`

— Method`simpson(solution; threshold = 0, kwargs...)`

Computes the average Simpson diversity index, i.e. the second Hill number, over the `last`

timesteps.

`kwargs...`

arguments are forwarded to `extract_last_timesteps`

. See `extract_last_timesteps`

for the argument details.

Can also handle a vector, e.g. simpson(n::AbstractVector; threshold = 0)

**Reference**

https://en.wikipedia.org/wiki/Diversity*index#Simpson*index

`EcologicalNetworksDynamics.Internals.simulate`

— Method```
simulate(
params::ModelParameters,
B0::AbstractVector;
N0::AbstractVector = nothing,
alg = nothing,
t0::Number = 0,
tmax::Number = 500,
extinction_threshold::Union{Number,AbstractVector} = 1e-5,
verbose = true,
callback = CallbackSet(
TerminateSteadyState(1e-6, 1e-4),
ExtinctionCallback(extinction_threshold, length(B0), verbose),
),
diff_code_data = (dBdt!, params),
kwargs...,
)
```

Run biomass dynamics simulation, given model parameters (`params`

) and the initial biomass (`B0`

). You can choose your solver algorithm by specifying to the `alg`

keyword.

The dynamic is solved between t=`t0`

and t=`tmax`

.

By default, we give the following callbacks to `solve()`

:

`TerminateSteadyState`

(from DiffEqCallbacks) which ends the simulation when a steady state is reached`ExtinctionCallback`

which extinguishes species whose biomass goes under the`extinction_threshold`

(either a number or a vector, see`ExtinctionCallback`

).

You are free to provide other callbacks, either by changing the parameter values of the callbacks above, choosing other callbacks from DiffEqCallbacks or by creating you own callbacks.

Moreover, we use the `isoutofdomain`

argument of `solve()`

to reject time steps that lead to negative biomass values due to numerical error.

Extra performance may be achieved by providing specialized julia code to the `diff_code_data`

argument instead of the default, generic `EcologicalNetworksDynamics.dBdt!`

code. No need to write it yourself as `generate_dbdt`

does it for you.

Thanks to the extra keywords arguments `kwargs...`

, you have a direct access to the interface of solve. Thus you can directly specify keyword arguments of solve(), for instance if you want to hint toward stiff solver you can give as an argument to simulate `alg_hints=[:stiff]`

.

The output of this function is the result of `DifferentialEquations.solve()`

, to learn how to handle this output see Solution Handling.

If you are not interested in the biomass trajectories, but want to find directly the biomass at steady state see `find_steady_state`

.

**Examples**

```
julia> foodweb = FoodWeb([0 0; 1 0]); # create the foodweb
br = BioRates(foodweb; d = 0); # set natural death rate to 0
params = ModelParameters(foodweb; biorates = br);
B0 = [0.5, 0.5]; # set initial biomass
solution = simulate(params, B0); # run simulation
@assert is_terminated(solution) # => a steady state has been found
@assert solution[begin] == B0 # initial biomass equals initial biomass
@assert isapprox(solution[end], [0.188, 0.219]; atol = 1e-2) # steady state biomass
julia> using DifferentialEquations;
solution_custom_alg = simulate(params, B0; alg = BS5());
@assert isapprox(solution_custom_alg[end], [0.188, 0.219]; atol = 1e-2)
julia> import Logging; # TODO: remove when warnings removed from `generate_dbdt`.
# generate specialized code (same simulation)
xpr, data = Logging.with_logger(Logging.NullLogger()) do
generate_dbdt(params, :raw)
end;
solution = simulate(params, B0; diff_code_data = (eval(xpr), data));
@assert is_terminated(solution) # the same result is obtained, more efficiently.
@assert solution[begin] == B0
@assert isapprox(solution[end], [0.188, 0.219]; atol = 1e-2) # steady state biomass
julia> # Same with alternate style.
xpr, data = Logging.with_logger(Logging.NullLogger()) do
generate_dbdt(params, :compact)
end;
solution = simulate(params, B0; diff_code_data = (eval(xpr), data));
@assert is_terminated(solution)
@assert solution[begin] == B0
@assert isapprox(solution[end], [0.188, 0.219]; atol = 1e-2) # steady state biomass
```

By default, the extinction callback throw a message when a species goes extinct.

```
julia> foodweb = FoodWeb([0 0; 1 0]);
params = ModelParameters(foodweb; biorates = BioRates(foodweb; d = 0));
simulate(params, [0.5, 1e-12]; verbose = true); # default: a message is thrown
[ Info: Species [2] is exinct. t=0.12316364776188903
julia> simulate(params, [0.5, 1e-12]; verbose = true); # no message thrown
```

`EcologicalNetworksDynamics.Internals.species_cv`

— Method`species_cv(solution::Solution; threshold = 0, last = "10%", kwargs...)`

Computes the temporal coefficient of variation of species biomass and its average.

See `coefficient_of_variation`

for the details.

`EcologicalNetworksDynamics.Internals.species_indices`

— MethodSpecies indices of the given `network`

.

`EcologicalNetworksDynamics.Internals.species_persistence`

— Method`species_persistence(solution; kwargs...)`

Returns the average proportion of species having a biomass superior or equal to threshold over the `last`

timesteps.

`kwargs...`

arguments are forwarded to `extract_last_timesteps`

. See `extract_last_timesteps`

for the argument details.

When applied to a vector of biomass, e.g. `species_persistence(n::Vector; threshold = 0)`

, it returns the proportion of species which biomass is above `threshold`

.

**Examples**

```
julia> foodweb = FoodWeb([0 0; 0 0]; quiet = true);
params = ModelParameters(foodweb);
sim_two = simulate(params, [0.5, 0.5]);
species_persistence(sim_two; last = 1)
1.0
julia> sim_one = simulate(params, [0, 0.5]);
species_persistence(sim_one; last = 1)
0.5
julia> sim_zero = simulate(params, [0, 0]);
species_persistence(sim_zero; last = 1)
0.0
julia> species_persistence([0, 1])
0.5
julia> species_persistence([1, 1])
1.0
```

`EcologicalNetworksDynamics.Internals.symmetrize`

— Method`symmetrize(V)`

Add symmetric tuples from an asymmetric vector of tuples. A vector `V`

of tuples is said to be asymmetric ⟺ ((i,j) ∈ `V`

⇒ (j,i) ∉ `V`

).

See also `asymmetrize`

.

`EcologicalNetworksDynamics.Internals.synchrony`

— Method```
synchrony(
solution::Solution;
threshold = 0,
last = "10%",
corrected = true,
kwargs...,
)
```

Compute the synchrony of species biomass fluctuations following Loreau & de Mazancourt (2008).

See `coefficient_of_variation`

for the argument details.

**Reference**

Loreau, M., & de Mazancourt, C. (2008). Species Synchrony and Its Drivers : Neutral and Nonneutral Community Dynamics in Fluctuating Environments. The American Naturalist, 172(2), E48‑E66. https://doi.org/10.1086/589746

`EcologicalNetworksDynamics.Internals.top_predators`

— Method`top_predators(net::EcologicalNetwork)`

Top predator indices (species eaten by nobody) of the given `net`

work.

```
julia> foodweb = FoodWeb([1 => 2, 2 => 3])
top_predators(foodweb) == [1]
true
```

See also `trophic_levels`

and `trophic_classes`

.

`EcologicalNetworksDynamics.Internals.total_richness`

— Method`total_richness(p::ModelParameters)`

Total richness of the system defined as the sum of the number of species and the number of nutrients. If the model `p`

does not include nutrient dynamics for producer growth, then `total_richness`

is equivalent the species `richness`

.

```
julia> foodweb = FoodWeb([1 => 2]) # 1 eats 2.
producer_growth = LogisticGrowth(foodweb)
model = ModelParameters(foodweb; producer_growth)
total_richness(model) == richness(model) == 2
true
julia> producer_growth = NutrientIntake(foodweb; n_nutrients = 3) # Include nutrients.
model_with_nutrients = ModelParameters(foodweb; producer_growth)
S = richness(model_with_nutrients)
N = nutrient_richness(model_with_nutrients)
total_richness(model_with_nutrients) == S + N == 5 # 2 + 3
true
```

See also `richness`

and `nutrient_richness`

.

`EcologicalNetworksDynamics.Internals.trophic_classes`

— Method`trophic_classes(foodweb::FoodWeb)`

Categorize each species into three classes: either `producers`

, `intermediate_consumers`

or `top_predators`

.

```
julia> foodweb = FoodWeb([1 => 2, 2 => 3])
trophic_classes(foodweb)
(producers = [3], intermediate_consumers = [2], top_predators = [1])
```

To access the name of the trophic classes you can simply call `trophic_classes`

without any argument.

```
julia> trophic_classes()
(:producers, :intermediate_consumers, :top_predators)
```

See also `trophic_levels`

and `top_predators`

.

`EcologicalNetworksDynamics.Internals.trophic_levels`

— Method`trophic_levels(net::EcologicalNetwork)`

Trophic level of each species.

```
julia> foodweb = FoodWeb([1 => 2])
trophic_levels(foodweb) == [2.0, 1.0]
true
```

See also `top_predators`

and `trophic_classes`

.

`EcologicalNetworksDynamics.Internals.trophic_structure`

— Method`trophic_structure(solution; threshold = 0, idxs = nothing, kwargs...)`

Returns the maximum, mean and weighted mean trophic level averaged over the `last`

timesteps. It also returns the adjacency matrix containing only the living species and the vector of the living species at the last timestep.

`extract_last_timesteps`

. See `extract_last_timesteps`

for the argument details.

**Examples**

```
julia> foodweb = FoodWeb([0 0 0; 0 0 0; 1 1 0]);
params = ModelParameters(foodweb);
B0 = [0.5, 0.5, 0.5];
sol = simulate(params, B0; verbose = true);
three_sp = trophic_structure(sol; last = 10);
three_sp[(:max, :mean)]
(max = 2.0, mean = 1.3333333333333335)
julia> B0 = [0.5, 0.5, 0];
sol = simulate(params, B0; verbose = true);
no_consumer = trophic_structure(sol; last = 10);
no_consumer[(:max, :mean)]
(max = 1.0, mean = 1.0)
julia> foodweb = FoodWeb([0 0 0; 0 1 0; 1 1 0]; quiet = true);
params = ModelParameters(foodweb);
B0 = [0.5, 0.5, 0.5];
sol = simulate(params, B0; verbose = false);
sum(trophic_structure(sol; last = 1).alive_A) - sum(foodweb.A)
-2
```

`EcologicalNetworksDynamics.Internals.vertebrates`

— MethodReturn indices of the vertebrates of the network (`net`

).

`EcologicalNetworksDynamics.Internals.weighted_mean_trophic_level`

— Method```
max_trophic_level(solution::Solution; threshold = 0, kwargs...)
mean_trophic_level(solution::Solution; threshold = 0, kwargs...)
weighted_mean_trophic_level(solution::Solution; threshold = 0, kwargs...)
```

Return the aggregated trophic level over the `last`

timesteps, either with `max`

or `mean`

aggregation, or by the mean trophic level weighted by species biomasses.

`extract_last_timesteps`

. See `extract_last_timesteps`

for the argument details.

These functions also handle biomass vectors associated with a network, as well as a vector of trophic levels (See examples).

**Examples**

```
julia> A = [0 0; 1 0];
max_trophic_level(A)
2.0
julia> mean_trophic_level(A)
1.5
julia> bm = [1, 1];
max_trophic_level(bm, A)
2.0
julia> mean_trophic_level(bm, A)
1.5
julia> weighted_mean_trophic_level(bm, A)
1.5
julia> bm = [0, 1];
max_trophic_level(bm, A)
1.0
julia> mean_trophic_level(bm, A)
1.0
julia> weighted_mean_trophic_level(bm, A)
1.0
julia> foodweb = FoodWeb(A; quiet = true);
params = ModelParameters(foodweb);
B0 = [0.5, 0.5];
sol = simulate(params, B0; verbose = false);
max_trophic_level(sol)
2.0
julia> mean_trophic_level(sol)
1.5
julia> w = weighted_mean_trophic_level(sol);
round(w; digits = 2)
1.54
```

`EcologicalNetworksDynamics.Internals.xp_sum`

— MethodRepeat the given expression into terms of a sum, successively replacing `indexes`

in `term`

by elements in (zipped) `lists`

.

```
julia> import EcologicalNetworksDynamics.Internals: xp_sum
julia> xp_sum([:i], [[1, 2, 3]], :(u^i)) # Three terms.
:(u ^ 1 + u ^ 2 + u ^ 3)
julia> xp_sum([:i], [[1]], :(u^i)) # Single term.
:(u ^ 1)
julia> xp_sum([:i], [[]], :(u^i)) # No terms.
0
julia> xp_sum([:i, :j], [[:a, :b, :c], [5, 8, 13]], :(j * i)) # Zipped indices.
:(5a + 8b + 13c)
```

`EcologicalNetworksDynamics.Foodweb`

— TypeThe Foodweb component, aka. "Trophic layer", adds a set of trophic links connecting species in the model. This is one of the most structuring components.

**Blueprint Creation from Raw Links.**

From an adjacency list:

```
julia> fw = Foodweb([:a => :b, :b => [:c, :d]])
blueprint for Foodweb with 2 trophic links:
A:
:a eats :b
:b eats :c and :d
julia> Model(fw) # Automatically brings a 'Species' component.
Model with 2 components:
- Species: 4 (:a, :b, :c, :d)
- Foodweb: 3 links
julia> Model(Foodweb([4 => [2, 1], 2 => 1])) # From species indices.
Model with 2 components:
- Species: 4 (:s1, :s2, :s3, :s4)
- Foodweb: 3 links
```

From a matrix:

```
julia> fw = Foodweb([0 0 1; 1 0 1; 0 0 0])
blueprint for Foodweb with 3 trophic links:
A: 3×3 SparseArrays.SparseMatrixCSC{Bool, Int64} with 3 stored entries:
⋅ ⋅ 1
1 ⋅ 1
⋅ ⋅
julia> Model(fw)
Model with 2 components:
- Species: 3 (:s1, :s2, :s3)
- Foodweb: 3 links
```

**Blueprint Creation from Random Models.**

**Cascade** model: specify the desired number of species `S`

and connectance `C`

.

```
julia> using Random
Random.seed!(12)
fw = Foodweb(:cascade; S = 5, C = 0.2)
blueprint for Foodweb with 5 trophic links:
A: 5×5 SparseArrays.SparseMatrixCSC{Bool, Int64} with 5 stored entries:
⋅ 1 ⋅ 1 1
⋅ ⋅ ⋅ ⋅ ⋅
⋅ ⋅ ⋅ 1 ⋅
⋅ ⋅ ⋅ ⋅ 1
⋅ ⋅ ⋅ ⋅ ⋅
```

Random foodwebs are drawn until the desired connectance is obtained, within a tolerance level defaulted to `tol_C = 0.1 * C`

, modifiable as a keyword argument.

**Niche** model: either specify the connectance `C`

or number of links `L`

.

```
julia> fw = Foodweb(:niche; S = 5, C = 0.2) # From connectance.
blueprint for Foodweb with 5 trophic links:
A: 5×5 SparseArrays.SparseMatrixCSC{Bool, Int64} with 5 stored entries:
⋅ ⋅ ⋅ ⋅ ⋅
⋅ ⋅ ⋅ ⋅ ⋅
1 1 ⋅ ⋅ ⋅
⋅ ⋅ 1 ⋅ ⋅
1 1 ⋅ ⋅ ⋅
julia> fw = Foodweb(:niche; S = 5, L = 4) # From number of links.
blueprint for Foodweb with 4 trophic links:
A: 5×5 SparseArrays.SparseMatrixCSC{Bool, Int64} with 4 stored entries:
⋅ ⋅ ⋅ ⋅ ⋅
⋅ ⋅ ⋅ ⋅ ⋅
⋅ ⋅ ⋅ ⋅ ⋅
1 ⋅ ⋅ ⋅ ⋅
1 1 1 ⋅ ⋅
```

The default tolerance levels for the niche model are `tol_C = 0.1 * C`

and `tol_L = 0.1 * L`

, modifiable as keyword arguments.

For either random model, the following keyword arguments can also be specified:

`reject_cycles = false`

(default): raise to forbid trophic cycles.`reject_if_disconnected = true`

(default): lower to allow disconnected trophic networks.`max_iterations = 10^5`

(default): give up if no satisfying network can be found after this number of random trials.

**Properties.**

A model `m`

with a `Foodweb`

has the following properties.

`m.A`

or`m.trophic_links`

: a view into the matrix of trophic links.`m.n_trophic_links`

: the number of trophic links in the model.`m.trophic_levels`

: calculate the trophic level of every species in the model.- Distinguishing between
`producers`

(species without outgoing trophic links) and`consumers`

(species with outgoing trophic links):`m.{producers,consumers}_mask`

: a boolean vector to select either kind of species.`m.n_{producers,consumers}`

: count number of species of either kind.`is_{producer,consumer}(m, i)`

: check whether species`i`

(name or index) is of either kind.`m.{producers,consumer}_indices`

: iterate over either species kind indices.`m.{producers,consumer}_{sparse,dense}_index`

: obtain a $species\_name \mapsto species\_index$ mapping:- the
`sparse`

index yields indices valid within the whole collection of species. - the
`dense`

index yields indices only valid within the restricted collection of species of either kind.

- the

- Distinguishing betwen
`preys`

(species with incoming trophic links) and`tops`

predators (species without incoming trophic links) works the same way. `m.producers_links`

: boolean matrix highlighting potential links between producers.`m.herbivorous_links`

: highlight only consumer-to-producer trophic links.`m.carnivorous_links`

: highlight only consumer-to-consumer trophic links.

```
julia> m = Model(Foodweb([:a => :b, :b => [:c, :d], :d => :e]));
julia> m.n_trophic_links
4
julia> m.A
5×5 EcologicalNetworksDynamics.TrophicLinks:
0 1 0 0 0
0 0 1 1 0
0 0 0 0 0
0 0 0 0 1
0 0 0 0 0
julia> m.trophic_levels
5-element EcologicalNetworksDynamics.TrophicLevels:
3.5
2.5
1.0
2.0
1.0
julia> m.producers_mask
5-element EcologicalNetworksDynamics.ProducersMask:
0
0
1
0
1
julia> m.preys_mask
5-element EcologicalNetworksDynamics.PreysMask:
0
1
1
1
1
julia> m.n_producers, m.n_consumers
(2, 3)
julia> m.n_tops, m.n_preys
(1, 4)
julia> is_top(m, 1), is_top(m, 2)
(true, false)
julia> collect(m.consumers_indices)
3-element Vector{Int64}:
1
2
4
julia> m.producers_sparse_index
OrderedCollections.OrderedDict{Symbol, Int64} with 2 entries:
:c => 3
:e => 5
julia> m.producers_dense_index
OrderedCollections.OrderedDict{Symbol, Int64} with 2 entries:
:c => 1
:e => 2
julia> m.producers_links
5×5 EcologicalNetworksDynamics.ProducersLinks:
0 0 0 0 0
0 0 0 0 0
0 0 1 0 1
0 0 0 0 0
0 0 1 0 1
julia> m.herbivorous_links
5×5 EcologicalNetworksDynamics.HerbivorousLinks:
0 0 0 0 0
0 0 1 0 0
0 0 0 0 0
0 0 0 0 1
0 0 0 0 0
julia> m.carnivorous_links
5×5 EcologicalNetworksDynamics.CarnivorousLinks:
0 1 0 0 0
0 0 0 1 0
0 0 0 0 0
0 0 0 0 0
0 0 0 0 0
```

`EcologicalNetworksDynamics.Model`

— TypeModel is the main object that we hand out to user which contains all the information about the underlying ecological model.

**Create a Model**

The most straightforward way to create a model is to use `default_model`

. This function only requires you to specify the trophic network.

```
fw = [1 => 2, 2 => 3]
model = default_model(fw)
```

This function will help you to create a model with ease, however it relies on default values for the parameters, which are not always suitable for your specific case, even though extracted from the literature.

To create a model with custom parameters, you can pass other arguments to `default_model`

.

`model = default_model(fw, BodyMass(; Z = 100))`

For instance, the above example creates a model with a body mass distribution with a predator-prey mass ratio of 100.

It is also possible to create a model manually by adding the components one by one. First, create an empty model:

`m = Model()`

Then add your components one by one. Note that you have to add the components in the right order, as some components depend on others. Moreover, some components are mandatory. Specifically, you need to provide a food web, species body masses, a functional response, metabolic rates and a producer growth function.

```
m = Model()
m += Foodweb([3 => 2, 2 => 1])
m += ClassicResponse(; h = 2, M = BodyMass([0.1, 2, 3]))
m += LogisticGrowth(; r = 1, K = 10)
m += Metabolism(:Miele2019)
m += Mortality(0)
```

**Read and write properties of the model**

First all properties contained in the model can be listed with:

`properties(m) # Where m is a Model.`

Then, the value of a property can be read with `get_<X>`

where `X`

is the name of the property. For instance, to read mortality rates:

`get_mortality(m) # Equivalent to: m.mortality.`

You can also re-write properties of the model using `set_<X>!`

. However, not all properties can be re-written, because some of them are derived from the others. For instance, many parameters are derived from species body masses, therefore changing body masses would make the model inconsistent. However, terminal properties can be re-written, as the species metabolic rate.

`EcologicalNetworksDynamics.Species`

— TypeThe Species component adds the most basic nodes compartment into the model: species. There is one node per species, and every species is given a unique name and index. The species ordering specified in this compartment is the reference species ordering.

```
julia> sp = Species(["hen", "fox", "snake"])
blueprint for Species:
names: 3-element Vector{Symbol}:
:hen
:fox
:snake
julia> m = Model(sp)
Model with 1 component:
- Species: 3 (:hen, :fox, :snake)
julia> Model(Species(5)) # Default names generated.
Model with 1 component:
- Species: 5 (:s1, :s2, :s3, :s4, :s5)
```

Typically, the species component is implicitly brought by other blueprints.

```
julia> Model(Foodweb([:a => :b]))
Model with 2 components:
- Species: 2 (:a, :b)
- Foodweb: 1 link
julia> Model(BodyMass([4, 5, 6]))
Model with 2 components:
- Species: 3 (:s1, :s2, :s3)
- Body masses: [4.0, 5.0, 6.0]
```

The species component makes the following properties available to a model `m`

:

`m.S`

or`m.richness`

or`m.species_richness`

or`m.n_species`

: number of species in the model.`m.species_names`

: list of species name in reference order.`m.species_index`

: get a $species\_name \mapsto species\_index$ mapping.

```
julia> (m.S, m.richness, m.species_richness, m.n_species) # All aliases for the same thing.
(3, 3, 3, 3)
julia> m.species_names
3-element EcologicalNetworksDynamics.SpeciesNames:
:hen
:fox
:snake
julia> m.species_index
OrderedCollections.OrderedDict{Symbol, Int64} with 3 entries:
:hen => 1
:fox => 2
:snake => 3
```

`EcologicalNetworksDynamics.default_model`

— Method```
default_model(
foodweb::Foodweb;
kwargs...
```

)

Generate a model from a food web with parameters set to default values.

Let's first illustrate the use of `default_model`

with a simple example.

```
foodweb = Foodweb([1 => 2])
default_model(foodweb)
```

In this example, *all* parameters are set to default values, however for your needs, you can override any of the default parameters. For instance, if you want to override the default metabolic rate, you can do it as follows:

```
my_x = [0.0, 1.2] # One value per species.
default_model(foodweb, Metabolism(my_x))
```

`EcologicalNetworksDynamics.persistence`

— Method`persistence(vec::AbstractVector; threshold = 0)`

Fraction of alive species given a biomass vector `vec`

. See `richness`

for details.

`EcologicalNetworksDynamics.persistence`

— Method`persistence(solution::AbstractODESolution; threshold = 0)`

Fraction of alive species at each timestep of the simulation. See `richness`

for details.

**Examples**

```
julia> S = 20 # Initial number of species.
foodweb = Foodweb(:niche; S = 20, C = 0.1)
m = default_model(foodweb)
B0 = rand(S)
sol = simulate(m, B0)
all(persistence(sol) .== richness(sol) / S)
true
```

`EcologicalNetworksDynamics.richness`

— Method`richness(vec::AbstractVector; threshold = 0)`

Return the number of alive species given a biomass vector `vec`

. By default, species are considered extinct if their biomass is 0. But, this `threshold`

can be changed using the corresponding keyword argument.

**Examples**

```
julia> foodweb = Foodweb([0 0; 1 0])
m = default_model(foodweb)
B0 = [0.5, 0.5]
sol = simulate(m, B0)
richness(sol[end]) # Richness at the end of the simulation.
2.0
```

`EcologicalNetworksDynamics.richness`

— Methodrichness(solution::Solution; threshold = 0)

Return the number of alive species at each timestep of the simulation. `solution`

is the output of `simulate`

. By default, species are considered extinct if their biomass is 0. But, this `threshold`

can be changed using the corresponding keyword argument.

**Examples**

Let's start with a simple example where the richness remains constant:

```
julia> foodweb = Foodweb([0 0; 1 0])
m = default_model(foodweb)
B0 = [0.5, 0.5]
sol = simulate(m, B0)
richness_trajectory = richness(sol)
all(richness_trajectory .== 2) # At each timestep, there are 2 alive species.
true
```

Now let's assume that the producer is extinct at the beginning of the simulation, while its consumer is not. We expect to observe a decrease in richness from 1 to 0 over time.

```
julia> B0 = [0, 0.5] # The producer is extinct at the beginning.
sol = simulate(m, B0)
richness_trajectory = richness(sol)
richness_trajectory[1] == 1 && richness_trajectory[end] == 0
true
```

`EcologicalNetworksDynamics.shannon_diversity`

— Method`shannon_diversity(vec::AbstractVector; threshold = 0)`

Shannon diversitty index given a biomass vector `vec Shannon diversity is a measure of species diversity based on the entropy. According to the Shannon index, for a same number of species, the more evenly the biomass is distributed among them, the higher the diversity.

**Example**

We consider a simple example with 3 species, but different shannon diversity.

```
julia> s1 = shannon_diversity([1, 1, 1])
s2 = shannon_diversity([1, 1, 0.1])
s3 = shannon_diversity([1, 1, 0.01])
s1 > s2 > s3
true
```

We observe as we decrease the biomass of the third species, the shannon diversity tends to 2, as we tend towards an effective two-species community.

`EcologicalNetworksDynamics.shannon_diversity`

— Method`shannon_diversity(solution::AbstractODESolution; threshold = 0)`

Shannon diversity index at each timestep of the simulation. `solution`

is the output of `simulate`

. Shannon diversity is a measure of species diversity based on the entropy. According to the Shannon index, for a same number of species, the more evenly the biomass is distributed among them, the higher the diversity.

**Example**

We start a simple simulation with even biomass distribution, therefore we expect the Shannon diversity to decrease over time as the biomass of the species diverge from each other.

```
julia> foodweb = Foodweb([0 0; 1 0])
m = default_model(foodweb)
B0 = [0.5, 0.5] # Even biomass, maximal shannon diversity.
sol = simulate(m, B0)
shannon_trajectory = shannon_diversity(sol)
biomass_trajectory[1] > biomass_trajectory[end]
true
```

`EcologicalNetworksDynamics.total_biomass`

— Method`total_biomass(vec::AbstractVector)`

Total biomass of a community given a biomass vector `vec`

.

**Examples**

```
julia> total_biomass([0.5, 1.5]) # 0.5 + 1.5 = 2.0
2.0
```

`EcologicalNetworksDynamics.total_biomass`

— Method`total_biomass(solution::AbstractODESolution)`

Total biomass of a community at each timestep of the simulation. `solution`

is the output of `simulate`

.

**Example**

Let's consider a consumer feeding on a producer, and let's start the simulation with the producer extinction so we can observe the consumer's biomass decrease over time.

```
julia> foodweb = Foodweb([0 0; 1 0])
m = default_model(foodweb)
B0 = [0, 0.5] # The producer is extinct at the beginning.
sol = simulate(m, B0)
biomass_trajectory = total_biomass(sol)
biomass_trajectory[1] == 0.5 && biomass_trajectory[end] == 0
true
```