Function to calculate trophic level

[gustavdelius]

I propose to add a function getTrophicLevel() that returns a matrix (species x size).

The trophic level of an individual is defined as 1 more than the average trophic level of all the prey items they have consumed during their life up to the current moment. The average is taken as a weighted average, weighted by the biomass of the prey. The trophic level of the primary resource is set to 0.

Traditionally one consider the simple case where all individuals of a species have the same diet composition $D_{ij}$ throughout their life, where $D_{ij}$ is the proportion of the total biomass intake of an individual of species $i$ that comes from prey from species $j$. Then we have that same trophic level for all individuals of species $i$, given by

$T_i = \sum_{j}D_{ij}T_j + 1.$

In vector notation this is just $\mathbf{T}=D\mathbf{T}+\mathbf{1}$ and thus $\mathbf{T}$ is obtained as the solution of

$(1-D)\mathbf{T}=\mathbf{1}.$

Thus in this case the trophic level is easy to calculate.

The definition of the trophic level of an individual becomes a bit more complicated if the diet changes as an individual grows up, as it does in mizer. We now need to integrate over the entire life of the predator until it reaches its current size.

We'll discuss this only in the simple case where the system is in a steady state so that the growth rate and the prey density depend on size only and not explicitly on time. We can replace the integral over time since birth by an integral over size since birth, using that

$dt = \frac{dw}{dw/dt}=\frac{dw}{g(w)}$

where $g(w)$ is the birth rate.
Let us denote by $r_{ij}(w,w_p)$ the rate density at which a predator of species $i$ and weight $w$ consumes biomass from prey species $j$ at weight $w_p$, which in our current scenario is time-independent. Then the trophic level of such a predator is

$T_i(w)=\frac{\int_{w_0}^w \frac{1}{g_i(w')}\sum_j \int r_{ij}(w',w_p)T_j(w_p) dw_p dw'}{\int_{w_0}^w \frac{1}{g_i(w')}\sum_j \int r_{ij}(w',w_p) dw_p dw'}+ 1.$

The denominator is the total biomass that the predator has consumed in their lifetime. Note that we need to include the resource as one of the prey "species". The resource will not contribute to the numerator because its trophic level is taken to be $0$ at all sizes, but it does contribute to the denominator.
We can again view the above equation as a vector equation $(1-D)\mathbf{T}=\mathbf{1}$ where now the entries of $\mathbf{T}$ are indexed by $(i,w)$ and

$D_{(i,w),(j,w_p)}=\frac{\int_{w_0}^w \frac{1}{g_i(w')} r_{ij}(w',w_p) dw'}{\int_{w_0}^w \frac{1}{g_i(w')}\sum_j \int r_{ij}(w',w_p) dw_p dw'}$

The matrix multiplication now involves both a sum over the prey species and an integral over the prey sizes. If we discretise the weight then $D$ becomes an $n\times n$ matrix where $n$ is the product of the number of species times the number of size grid points. The resulting $n$-dimensional linear system is still quite fast to solve.

In mizer the consumption rate density is given by

$r_{ij}(w,w_p)=\theta_{ij}\gamma_i(w)(1-f_i(w))\phi_i(w/w_p)N_j(w_p)w_p,$

where $\theta$ is the interaction matrix, $\gamma$ is the search volume, $f$ is the feeding level, $\phi$ is the predation kernel and $N$ is the abundance density.

One could go on and define the trophic level of a species as the average trophic level of all its individuals. This might be useful in order to compare with the trophic level of a species in a non-size-resolved model. One would probably want this to be a weighted average where the trophic level of each individual is weighted by the consumption rate of that individual:

$T_i = \frac{\int r_i(w)N_i(w)T_i(w)dw}{\int r_i(w)N_i(w)dw},$

where the consumption rate of an individual of species $i$ and weight $w$ is

$r_i(w)=\sum_j\int r_{ij}(w,w_p)dw_p.$

This could be returned by a function getTrophicLevelBySpecies().