Food Web Analysis

ZHAW Institute for Natural Resource Sciences
Module Applied Geoinformatics

Summary

As ecosystems are threatened by anthropogenic activities, food webs provide a way to comprehend these complex networks and potentially find ways to help sustain their function. Using methods from the analysis of geographic networks the keystone species of a food web have been determined and the stability of the network has been assessed. Six different keystone species were found, depending on the centrality index used. The food web has been visualised to give an overview of the network functions and to compare the indices differences.

KEYWORDS: Food web, Network, Keystone species, Centrality index, Connectivity

1 Introduction

1.1 Food webs

Food webs are a form of ecological networks that represent the food chains of an ecosystem. These networks include no geographic relation, however the analysis of their structure and properties using python shows great similarity to analysing geographic networks. Though they do not include coordinates, the structure is essentially the same. Like spatial networks, food webs consist of nodes and edges. Likewise, connectivity, centrality, and shortest paths are all important aspects for their analysis. Each node represents a species, while the edges represent interactions between said species. The edges can either be weighted for the species interaction magnitude, or be of simple binary nature. These models are used to better understand ecosystems, to assess their quality, or in further modelling approaches to detect possible changes in ecosystem structure, due to factors such as an introduced species or environmental factors. In so called “removal analyses” the effect of removing certain species can be modelled (Hui, 2012; Pimm et al., 1991; Solé & Montoya, 2001) Ecopath is a modelling software based on mass-balance and a set of mathematical formulas, which can be used to construct a food web. It generates a snapshot of the ecosystem structure, which is a static model. Further studies of the ecosystem use Ecosim, which includes the temporal changes for the model, and Ecospace, which includes the spatial heterogeneity of the model (Christensen et al., 2005).

1.2 Indices

To describe networks a variety of indices can be calculated. Some indices describe the network in its entirety, and others inspect individual nodes. Centrality indices belong to the latter sort and play a big role in determining keystone species. These are the species that are most crucial for sustaining network stability. The search for the most accurate indices for assessing keystone species has been subject to a variety of studies in the past decades (Solé & Montoya, 2001), (Jordán et al., 2007). While a consensus has not been reached, methods to assess the stability of indices have been developed (Fedor & Vasas, 2009). The indices approach keystone-ness of a species in a network in different ways and show different aspects of their influence on the ecosystem. Following formulas for indices have been recalled from networkx documentation (Hagberg et al., 2008), and web of life documentation (Vindigni, 2022).

1.2.1 Degree centrality

Degree centrality index describes the number of connections a node has in relation to the total possible connections. A species with a high degree centrality acts as a local hub. This means that a large number of species are directly affected by its abundance. What the index doesn’t indicate is whether the node has a large impact on the entire network, since only immediate connections are considered (Rocchi et al., 2017). Further, it does not consider weighted edges and treats every edge in a binary manner. Degree centrality can be calculated as shown in (1). \(d_v\) describes the degree of the node in question and \(N\) describes the number of total nodes (Rocchi et al., 2017).

\[\begin{align} c_d (v)= \frac{d_v}{N-1} \tag{1}\label{eq:1} \\ \end{align}\]

1.2.2 Betweenness centrality

Betweenness centrality index determines how often a node is included in the shortest paths between every node in the network. In a food web, nodes with high betweenness centrality indices act as transmitting paths to spread effects through the network (Rocchi et al., 2017). The calculation of betweenness centrality is defined as shown in (2).

\[\begin{align} c_b (v)= \sum_{s,t,\epsilon V}^{} \frac{\sigma (s,t\mid v)}{\sigma (s,t)} \tag{2}\label{eq:2} \\ \end{align}\]

1.2.3 Closeness centrality

A nodes closeness centrality is biggest when its shortest paths to all other nodes are of minimal length. For food webs that means that the effect caused by the species abundance is spread most rapidly across the whole network (Rocchi et al., 2017). For the calculation, the improved formula from Wasserman and Faust has been used (Wasserman and Faust as cited in Hagberg et al., 2008). The formula is shown in (3). \(n\) describes the number of reachable nodes from the node in question, \(d(v,u)\) describes the shortest path between nodes \(v\) and \(u\).

\[\begin{align} C_{WF}(v)=\frac{n-1}{N-1}\frac{n-1}{\sum_{v=1}^{n-1}d(v,u)} \tag{3}\label{eq:3} \\ \end{align}\]

1.2.4 In-degree

The in-degree is a measurement for directed networks. It assesses how many edges are pointed towards a particular node. In a food web, a high in-degree is characteristic for generalist predators that feed on many different species of prey (Rocchi et al., 2017).

1.2.5 Out-degree

Like the in-degree, the out-degree is used in directed networks. It measures how many edges are pointed away from the node in question. In a food web a high out-degree is characteristic for important prey species that serve as food source for many different predator species (Rocchi et al., 2017).

1.2.6 Trophic levels

The trophic level of every species can be calculated using the formula of (4). \(k_v^in\) describes the in-degree of a node, \(s_j\) describes the shortest path to a primary producer node and \(a_{vj}\) describes

\[\begin{align} s(v)=1+\frac{1}{k_{v}^{in}}\sum_{j}^{}a_{vj}s_j \tag{4}\label{eq:4} \\ \end{align}\]

1.2.7 Connectivity

The connectivity index is a measure describing the entire network and not an individual node. A highly connected network is more stable and can therefore better withstand disturbance. For food webs, the connectivity can be calculated as shown in (5). \(n_l\) describes the number of links and \(N\) describes the number of nodes.

\[\begin{align} C(G)=\frac{2n_l}{N(N-1)} \tag{5}\label{eq:5} \\ \end{align}\]

2 Methods

2.1 Data source

The food web from the Mondego Estuary in Portugal has been chosen for the analysis, taken from web-of-life, run by the University of Zurich (Vindigni, 2022). This dataset consists of 37 species and 242 interactions. It stems from the work of Baeta et al. who modelled the ecosystem of the estuary to assess its response to various external factors (Baeta et al., 2011).

2.2 Setup

To examine the data, python (version 3.10.6) has been used in combination with the jupyter notebook environment (version 6.4.12). The python library networkx was used to create the network graph and calculate its indices. The python library bokeh was used for visualization.

2.3 Analysis

The data was available in the form of a weighted neighborhood matrix. To construct a graph with networkx, a list of nodes and a list of edges had to be provided. Those were extracted from the matrix by accessing the respective columns and using a for loop. First an empty graph model in the form of a “MultiDiGraph” (directed graph with self loops and parallel edges) was created. Self loops had to be allowed because in some cases, species may feed on its own kind. Parallel edges can also be possible in food webs, because two species may feed on each other. Then, the nodes and edges were added, also providing the weight of the edges in the same step. For calculating indices, networkx provided all the functions required. The indices were saved in variables and later accessed when visualising. Bokeh does not accept graphs unless the node names are given as integer types. Therefore, the species names were changed to numbers.

3 Results

The benthic ecosystem of the Mondego estuary consists of three trophic levels. It shows a connectivity of 0.354. Five different indices have been calculated in the search for keystone species. The values of these calculations can be seen in Table 1. A visualisation of the network is given in Figure 7.

3.1 Import Data

At first the weighted neighborhood matrix (Table 1) has been imported.

Code

import pandas as pd
foodweb = pd.read_csv("FW_016_01.csv")
foodweb

Table 1: weighted neighborhood matrix

	Unnamed: 0	Alkmaria romijni	Ampithoe valida	Aonides oxycephala	Capitella capitata	Carcinus maenas	Cerastoderma edule	Chaetozone setosa	Crangon crangon	Cyathura carinata	...	Lumbrineris impatiens	Melita palmata	Mytilus sp1 FW_016_01	Nemertini	Nephthys sp1 FW_016_01	Oligochaeta	Pygospio elegans	Scrobicularia plana	Streblospio shrubsolii	Zooplankton
0	Alkmaria romijni	0.0	0.0	0	0	0.000	0.0	0	0.003	0.020	...	0.000	0.0	0.0	0.000	0.001	0	0.0	0.0	0	0.00
1	Ampithoe valida	0.0	0.0	0	0	0.040	0.0	0	0.037	0.050	...	0.001	0.0	0.0	0.050	0.006	0	0.0	0.0	0	0.00
2	Aonides oxycephala	0.0	0.0	0	0	0.001	0.0	0	0.001	0.001	...	0.001	0.0	0.0	0.005	0.001	0	0.0	0.0	0	0.00
3	Capitella capitata	0.0	0.0	0	0	0.000	0.0	0	0.000	0.004	...	0.002	0.0	0.0	0.000	0.012	0	0.0	0.0	0	0.00
4	Carcinus maenas	0.0	0.0	0	0	0.060	0.0	0	0.010	0.055	...	0.000	0.0	0.0	0.000	0.000	0	0.0	0.0	0	0.00
5	Cerastoderma edule	0.0	0.0	0	0	0.020	0.0	0	0.030	0.000	...	0.034	0.0	0.0	0.040	0.050	0	0.0	0.0	0	0.00
6	Chaetozone setosa	0.0	0.0	0	0	0.018	0.0	0	0.013	0.000	...	0.001	0.0	0.0	0.030	0.010	0	0.0	0.0	0	0.00
7	Crangon crangon	0.0	0.0	0	0	0.140	0.0	0	0.010	0.000	...	0.000	0.0	0.0	0.000	0.000	0	0.0	0.0	0	0.00
8	Cyathura carinata	0.0	0.0	0	0	0.030	0.0	0	0.040	0.000	...	0.010	0.0	0.0	0.030	0.000	0	0.0	0.0	0	0.00
9	Detritus	0.9	0.7	1	1	0.150	0.1	1	0.300	0.700	...	0.600	0.7	0.1	0.200	0.550	1	0.4	0.6	1	0.00
10	Diopatra neapolitana	0.0	0.0	0	0	0.001	0.0	0	0.001	0.000	...	0.000	0.0	0.0	0.005	0.001	0	0.0	0.0	0	0.00
11	Eteone flava	0.0	0.0	0	0	0.000	0.0	0	0.001	0.000	...	0.000	0.0	0.0	0.001	0.000	0	0.0	0.0	0	0.00
12	Gibulla umbilicalis	0.0	0.0	0	0	0.001	0.0	0	0.001	0.000	...	0.020	0.0	0.0	0.005	0.000	0	0.0	0.0	0	0.00
13	Glycera tridactyla	0.0	0.0	0	0	0.020	0.0	0	0.010	0.009	...	0.002	0.0	0.0	0.040	0.024	0	0.0	0.0	0	0.00
14	Green macroalgae	0.0	0.2	0	0	0.040	0.0	0	0.050	0.000	...	0.000	0.2	0.0	0.000	0.000	0	0.0	0.0	0	0.00
15	Haminoea hydatis	0.0	0.0	0	0	0.001	0.0	0	0.000	0.000	...	0.014	0.0	0.0	0.005	0.000	0	0.0	0.0	0	0.00
16	Hediste diversicolor	0.0	0.0	0	0	0.080	0.0	0	0.090	0.050	...	0.004	0.0	0.0	0.230	0.080	0	0.0	0.0	0	0.00
17	Heteromastus filiformis	0.0	0.0	0	0	0.005	0.0	0	0.040	0.030	...	0.005	0.0	0.0	0.077	0.100	0	0.0	0.0	0	0.00
18	Hydrobia ulvae	0.0	0.0	0	0	0.030	0.0	0	0.070	0.000	...	0.070	0.0	0.0	0.050	0.000	0	0.0	0.0	0	0.00
19	Idotea chelipes	0.0	0.0	0	0	0.005	0.0	0	0.005	0.010	...	0.001	0.0	0.0	0.010	0.000	0	0.0	0.0	0	0.00
20	Imports	0.0	0.0	0	0	0.180	0.0	0	0.060	0.000	...	0.000	0.0	0.0	0.000	0.000	0	0.0	0.0	0	0.00
21	Lagis koreni	0.0	0.0	0	0	0.001	0.0	0	0.002	0.000	...	0.000	0.0	0.0	0.005	0.000	0	0.0	0.0	0	0.00
22	Lekanesphaera levii	0.0	0.0	0	0	0.000	0.0	0	0.000	0.001	...	0.000	0.0	0.0	0.001	0.000	0	0.0	0.0	0	0.00
23	Littorina littorea	0.0	0.0	0	0	0.010	0.0	0	0.010	0.000	...	0.035	0.0	0.0	0.010	0.000	0	0.0	0.0	0	0.00
24	Lumbrineris impatiens	0.0	0.0	0	0	0.010	0.0	0	0.015	0.000	...	0.000	0.0	0.0	0.020	0.002	0	0.0	0.0	0	0.00
25	Melita palmata	0.0	0.0	0	0	0.030	0.0	0	0.030	0.045	...	0.002	0.0	0.0	0.040	0.008	0	0.0	0.0	0	0.00
26	Microphytobenthos	0.0	0.1	0	0	0.000	0.1	0	0.000	0.000	...	0.050	0.1	0.1	0.000	0.050	0	0.0	0.2	0	0.00
27	Mytilus sp1 FW_016_01	0.0	0.0	0	0	0.001	0.0	0	0.005	0.000	...	0.000	0.0	0.0	0.000	0.000	0	0.0	0.0	0	0.00
28	Nemertini	0.0	0.0	0	0	0.000	0.0	0	0.000	0.000	...	0.001	0.0	0.0	0.050	0.000	0	0.0	0.0	0	0.00
29	Nephthys sp1 FW_016_01	0.0	0.0	0	0	0.005	0.0	0	0.001	0.000	...	0.001	0.0	0.0	0.005	0.001	0	0.0	0.0	0	0.00
30	Oligochaeta	0.0	0.0	0	0	0.001	0.0	0	0.005	0.020	...	0.010	0.0	0.0	0.050	0.080	0	0.0	0.0	0	0.00
31	Phytoplankton	0.1	0.0	0	0	0.000	0.8	0	0.000	0.000	...	0.000	0.0	0.8	0.000	0.000	0	0.6	0.2	0	0.98
32	Pygospio elegans	0.0	0.0	0	0	0.000	0.0	0	0.000	0.000	...	0.000	0.0	0.0	0.001	0.000	0	0.0	0.0	0	0.00
33	Scrobicularia plana	0.0	0.0	0	0	0.060	0.0	0	0.060	0.000	...	0.035	0.0	0.0	0.030	0.024	0	0.0	0.0	0	0.00
34	Streblospio shrubsolii	0.0	0.0	0	0	0.000	0.0	0	0.000	0.005	...	0.001	0.0	0.0	0.005	0.000	0	0.0	0.0	0	0.00
35	Zooplankton	0.0	0.0	0	0	0.060	0.0	0	0.100	0.000	...	0.000	0.0	0.0	0.005	0.000	0	0.0	0.0	0	0.02
36	Zostera noltii	0.0	0.0	0	0	0.000	0.0	0	0.000	0.000	...	0.100	0.0	0.0	0.000	0.000	0	0.0	0.0	0	0.00

37 rows × 32 columns

3.2 Node list

Code

nodes_wduplicates = list(foodweb.iloc[:,0]) + list(foodweb.iloc[:, 1:-3].columns)
nodes = list(dict.fromkeys(nodes_wduplicates))

3.3 Edge list

For every edge in the network the connected nodes and weight of the link had to be extracted to (Table 2).

Code

u = list(range(0,len(foodweb)))
v = list(range(1,len(foodweb.columns)-3))

column_link = []
row_link = []
strength_link = []

for i in u:
    for j in v:
        strength_val = foodweb.iloc[i,j]
        if strength_val > 0:
            column_link.append(j)
            row_link.append(i)
            strength_link.append(strength_val)
            
column_link_names = list(foodweb.columns[column_link])
row_link_names = list(foodweb.iloc[row_link,0])

links = pd.DataFrame(list(zip(row_link_names,column_link_names,strength_link)))
links

Table 2: edge list with weight

	0	1	2
0	Alkmaria romijni	Crangon crangon	0.003
1	Alkmaria romijni	Cyathura carinata	0.020
2	Alkmaria romijni	Glycera tridactyla	0.010
3	Alkmaria romijni	Nephthys sp1 FW_016_01	0.001
4	Ampithoe valida	Carcinus maenas	0.040
...	...	...	...
231	Zooplankton	Crangon crangon	0.100
232	Zooplankton	Nemertini	0.005
233	Zostera noltii	Idotea chelipes	0.050
234	Zostera noltii	Littorina littorea	0.200
235	Zostera noltii	Lumbrineris impatiens	0.100

236 rows × 3 columns

3.4 Network graph

The graph has been constructed using networkx. In (Figure 1) a first visualisation of the network can be seen.

Code

import networkx as nx

G = nx.MultiDiGraph()

for node in nodes:
    G.add_node(node, name = node)

for link in list(links.index):
    G.add_edge(links.iloc[link,0], links.iloc[link,1], weight = links.iloc[link,2])

nx.draw(G, with_labels = True)

Figure 1: networkx graph

3.5 Converting to integer node names

Code

Gint = nx.convert_node_labels_to_integers(G, ordering="default")

3.6 Calculating indices

The indices have been calculated using networkx and summarised in (Table 3).

Code

# degree centrality
degree = nx.degree_centrality(Gint)
degree_list = list(degree.values())
indegree = nx.in_degree_centrality(Gint)
indegree_list = list(indegree.values())
outdegree = nx.out_degree_centrality(Gint)
outdegree_list = list(outdegree.values())

# closeness centrality
closeness = nx.closeness_centrality(Gint,distance='weight',wf_improved=True)
closeness_list = list(closeness.values())

# betweenness centrality
betweenness = nx.betweenness_centrality(Gint, weight='weight')
betweenness_list = list(betweenness.values())

#trophic levels
trophic = nx.trophic_levels(Gint,weight='weight')

indices = pd.DataFrame({'dc':pd.Series(degree), 'cc':pd.Series(closeness), 'bc':pd.Series(betweenness), 'in-degree':pd.Series(indegree), 'out-degree':pd.Series(outdegree), 'trophic':pd.Series(trophic)})
indices

C:\Users\LukasMeyer\AppData\Local\Programs\ArcGIS\Pro\bin\Python\envs\agi-env-2\lib\site-packages\networkx\algorithms\centrality\trophic.py:45: FutureWarning: adjacency_matrix will return a scipy.sparse array instead of a matrix in Networkx 3.0.
  a = nx.adjacency_matrix(G, weight=weight).T.toarray()

Table 3: indices for nodes

	dc	cc	bc	in-degree	out-degree	trophic
0	0.166667	0.111111	0.000000	0.055556	0.111111	2.000000
1	0.333333	0.250000	0.000661	0.083333	0.250000	2.000000
2	0.277778	0.027778	0.000000	0.027778	0.250000	2.000000
3	0.166667	0.027778	0.000000	0.027778	0.138889	2.000000
4	0.833333	42.056075	0.017857	0.750000	0.083333	2.677543
5	0.305556	0.250000	0.001852	0.083333	0.222222	2.000000
6	0.250000	0.027778	0.000000	0.027778	0.222222	2.000000
7	0.833333	63.380282	0.000000	0.777778	0.055556	2.498839
8	0.500000	27.777778	0.019048	0.388889	0.111111	2.357648
9	0.777778	0.000000	0.000000	0.000000	0.777778	1.000000
10	0.416667	3.260031	0.011905	0.277778	0.138889	2.180000
11	0.277778	12.491325	0.005556	0.194444	0.083333	2.784911
12	0.194444	0.250000	0.002778	0.083333	0.111111	2.000000
13	0.888889	31.746032	0.047619	0.638889	0.250000	2.707078
14	0.222222	0.000000	0.000000	0.000000	0.222222	1.000000
15	0.166667	0.250000	0.002778	0.083333	0.083333	2.000000
16	0.722222	32.967033	0.007143	0.500000	0.222222	2.330379
17	0.305556	0.027778	0.000000	0.027778	0.277778	2.000000
18	0.250000	0.111111	0.000000	0.055556	0.194444	2.000000
19	0.333333	1.000000	0.007143	0.166667	0.166667	2.250000
20	0.055556	0.000000	0.000000	0.000000	0.055556	1.000000
21	0.194444	0.027778	0.000000	0.027778	0.166667	2.000000
22	0.111111	0.111111	0.000000	0.055556	0.055556	2.000000
23	0.250000	0.444444	0.000000	0.111111	0.138889	2.000000
24	0.777778	34.220532	0.096429	0.638889	0.138889	2.221915
25	0.333333	0.250000	0.000661	0.083333	0.250000	2.000000
26	0.361111	0.000000	0.000000	0.000000	0.361111	1.000000
27	0.138889	0.250000	0.000794	0.083333	0.055556	2.000000
28	0.861111	30.176027	0.045238	0.750000	0.111111	2.932342
29	0.666667	32.490975	0.094841	0.472222	0.194444	2.420444
30	0.277778	0.027778	0.000000	0.027778	0.250000	2.000000
31	0.138889	0.000000	0.000000	0.000000	0.138889	1.000000
32	0.138889	0.111111	0.000000	0.055556	0.083333	2.000000
33	0.222222	0.000000	0.000000	0.000000	0.222222	1.000000
34	0.194444	0.000000	0.000000	0.000000	0.194444	1.000000
35	0.083333	0.000000	0.000000	0.000000	0.083333	1.000000
36	0.083333	0.000000	0.000000	0.000000	0.083333	1.000000

3.7 Ranking

The species have been ranked for every index calculated as it can be seen in (Table 4).

Code

keystone_rank = pd.DataFrame()
keystone_rank['dc_rank'] = indices['dc'].rank(ascending=False)
keystone_rank['cc_rank'] = indices['cc'].rank(ascending=False)
keystone_rank['bc_rank'] = indices['bc'].rank(ascending=False)
keystone_rank['indegree_rank'] = indices['in-degree'].rank(ascending=False)
keystone_rank['outdegree_rank'] = indices['out-degree'].rank(ascending=False)
keystone_rank['trophic_rank'] = indices['trophic'].rank(ascending=False)
keystone_rank

Table 4: nodes ranked by indices

	dc_rank	cc_rank	bc_rank	indegree_rank	outdegree_rank	trophic_rank
0	29.0	20.5	27.0	20.5	25.5	20.0
1	13.0	15.5	15.5	15.5	6.0	20.0
2	18.0	25.5	27.0	25.5	6.0	20.0
3	29.0	25.5	27.0	25.5	21.0	20.0
4	3.5	2.0	6.0	2.5	30.5	4.0
5	15.5	15.5	13.0	15.5	11.0	20.0
6	21.0	25.5	27.0	25.5	11.0	20.0
7	3.5	1.0	27.0	1.0	35.5	5.0
8	9.0	8.0	5.0	8.0	25.5	7.0
9	5.5	33.0	27.0	33.0	1.0	33.0
10	10.0	10.0	7.0	9.0	21.0	11.0
11	18.0	9.0	10.0	10.0	30.5	2.0
12	26.0	15.5	11.5	15.5	25.5	20.0
13	1.0	6.0	3.0	4.5	6.0	3.0
14	23.5	33.0	27.0	33.0	11.0	33.0
15	29.0	15.5	11.5	15.5	30.5	20.0
16	7.0	4.0	8.5	6.0	11.0	8.0
17	15.5	25.5	27.0	25.5	3.0	20.0
18	21.0	20.5	27.0	20.5	15.0	20.0
19	13.0	11.0	8.5	11.0	17.5	9.0
20	37.0	33.0	27.0	33.0	35.5	33.0
21	26.0	25.5	27.0	25.5	17.5	20.0
22	34.0	20.5	27.0	20.5	35.5	20.0
23	21.0	12.0	27.0	12.0	21.0	20.0
24	5.5	3.0	1.0	4.5	21.0	10.0
25	13.0	15.5	15.5	15.5	6.0	20.0
26	11.0	33.0	27.0	33.0	2.0	33.0
27	32.0	15.5	14.0	15.5	35.5	20.0
28	2.0	7.0	4.0	2.5	25.5	1.0
29	8.0	5.0	2.0	7.0	15.0	6.0
30	18.0	25.5	27.0	25.5	6.0	20.0
31	32.0	33.0	27.0	33.0	21.0	33.0
32	32.0	20.5	27.0	20.5	30.5	20.0
33	23.5	33.0	27.0	33.0	11.0	33.0
34	26.0	33.0	27.0	33.0	15.0	33.0
35	35.5	33.0	27.0	33.0	30.5	33.0
36	35.5	33.0	27.0	33.0	30.5	33.0

3.8 Connectivity

Code

connectivity = (2* len(Gint.edges())) / (len(nodes) * (len(nodes)-1))
connectivity

0.35435435435435436

3.9 Visualising graph

The network has been visualised using different indices for depicting the nodes. For (Figure 2), (Figure 3) and (Figure 4) the size and colour of the nodes represent the respective index value.

Code

from bokeh.models import Circle, MultiLine, StaticLayoutProvider
from bokeh.plotting import figure, from_networkx, show
from bokeh.io import output_notebook
from bokeh.transform import linear_cmap
from bokeh.palettes import Plasma8

p = figure(width=650, height=650, x_range=(-3, 3), y_range=(-3, 3),
           tools="hover,wheel_zoom,pan,reset", toolbar_location='above',
           x_axis_location=None, y_axis_location=None,
           title=None)

p.grid.grid_line_color = None

graph = from_networkx(Gint, nx.spring_layout, scale=2, center=(0,0))

graph.node_renderer.data_source.add([i for i in degree_list], 'degree')
graph.node_renderer.data_source.add([i*25+5 for i in degree_list], 'size')
graph.node_renderer.data_source.add([i for i in closeness_list], 'closeness')
graph.node_renderer.data_source.add([i for i in betweenness_list], 'betweenness')
graph.node_renderer.data_source.add([i for i in indegree_list], 'indegree')
graph.node_renderer.data_source.add([i for i in outdegree_list], 'outdegree')
graph.node_renderer.data_source.add([i for i in nodes], 'name')
graph.edge_renderer.data_source.add([i*3+1 for i in strength_link], 'weight')

graph.node_renderer.glyph = Circle(
    size='size', 
    fill_color=linear_cmap('degree', 'Plasma8', low=0, high=1))

graph.edge_renderer.glyph = MultiLine(line_width = 'weight')

p.hover.tooltips = [("index","@index"), ("name","@name"), ("degree","@degree"), ("closeness","@closeness"), ("betweenness","@betweenness"),
                    ("indegree","@indegree"), ("outdegree", "@outdegree")]

p.renderers.append(graph)

output_notebook()
show(p)

Loading BokehJS ...

Figure 2: Food web by degree centrality

Code

p = figure(width=650, height=650, x_range=(-3, 3), y_range=(-3, 3),
           tools="hover,wheel_zoom,pan,reset", toolbar_location='above',
           x_axis_location=None, y_axis_location=None,
           title=None)

p.grid.grid_line_color = None

graph = from_networkx(Gint, nx.spring_layout, scale=2, center=(0,0))

graph.node_renderer.data_source.add([i for i in degree_list], 'degree')
graph.node_renderer.data_source.add([i/4+10 for i in closeness_list], 'size')
graph.node_renderer.data_source.add([i for i in closeness_list], 'closeness')
graph.node_renderer.data_source.add([i for i in betweenness_list], 'betweenness')
graph.node_renderer.data_source.add([i for i in indegree_list], 'indegree')
graph.node_renderer.data_source.add([i for i in outdegree_list], 'outdegree')
graph.node_renderer.data_source.add([i for i in nodes], 'name')
graph.edge_renderer.data_source.add([i*3+1 for i in strength_link], 'weight')

graph.node_renderer.glyph = Circle(
    size='size', 
    fill_color=linear_cmap('closeness', 'Plasma8', low=0, high=1))

graph.edge_renderer.glyph = MultiLine(line_width = 'weight')

p.hover.tooltips = [("index","@index"), ("name","@name"), ("degree","@degree"), ("closeness","@closeness"), ("betweenness","@betweenness"),
                    ("indegree","@indegree"), ("outdegree", "@outdegree")]

p.renderers.append(graph)

output_notebook()
show(p)

Loading BokehJS ...

Figure 3: Food web by closeness centrality

Code

p = figure(width=650, height=650, x_range=(-3, 3), y_range=(-3, 3),
           tools="hover,wheel_zoom,pan,reset", toolbar_location='above',
           x_axis_location=None, y_axis_location=None,
           title=None)

p.grid.grid_line_color = None

graph = from_networkx(Gint, nx.spring_layout, scale=2, center=(0,0))

graph.node_renderer.data_source.add([i for i in degree_list], 'degree')
graph.node_renderer.data_source.add([i*200+10 for i in betweenness_list], 'size')
graph.node_renderer.data_source.add([i for i in closeness_list], 'closeness')
graph.node_renderer.data_source.add([i for i in betweenness_list], 'betweenness')
graph.node_renderer.data_source.add([i for i in indegree_list], 'indegree')
graph.node_renderer.data_source.add([i for i in outdegree_list], 'outdegree')
graph.node_renderer.data_source.add([i for i in nodes], 'name')
graph.edge_renderer.data_source.add([i*3+1 for i in strength_link], 'weight')

graph.node_renderer.glyph = Circle(
    size='size', 
    fill_color=linear_cmap('betweenness', 'Plasma8', low=0, high=0.1))

graph.edge_renderer.glyph = MultiLine(line_width = 'weight')

p.hover.tooltips = [("index","@index"), ("name","@name"), ("degree","@degree"), ("closeness","@closeness"), ("betweenness","@betweenness"),
                    ("indegree","@indegree"), ("outdegree", "@outdegree")]

p.renderers.append(graph)

output_notebook()
show(p)

Loading BokehJS ...

Figure 4: Food web by Betweenness centrality

4 Discussion

While every index provides a unique assessment of the species’ importance in the network, certain species are clearly more crucial for the ecosystem’s stability than others. Carcinus maenas (4), Crangon crangon (7), Glycera tridactyla (13), Lumbrineris impatiens (24), Nemertini (28) and Nephthys sp1 (29) could all be considered keystone species based on the different indices. In ecosystem research it has been suggested to shift from considering only one species as keystone species, to identifying a complex of keystone species as it proposes a more realistic and holistic approach (Ortiz et al., 2013). While analysing the network most keystone species are also observed to be categorised on a high trophic level. While it is true that apex predators play a vital role in ecosystems (Wallach et al., 2017), it should be noted that even primary producers can occupy keystone roles (Ortiz et al., 2013). Nemertini (28) shows the highest in-degree value and can thus be labeled as the most broadly feeding predator. However, because it is only determined to the genus level, this most likely comes to its advantage. Interestingly, Detritus (9) shows one of the highest degree centrality values in the ecosystem but closeness centrality and betweenness centrality show a value of zero. The same effect can be found for Microphytobenthos (26) and Phytoplankton (31). The problem is consistent for all primary producers in the ecosystem and arises from the directed nature of the network. Primary producers cannot own ingoing edges, so there is no way they are included in any shortest paths. Because degree centrality cannot interpret directed relationships it gives a value that is not zero.

Another interesting index to calculate would have been the keystone index. It uses both a top-down and bottom-up approach for calculating the species importance (Jordán as cited in Fedor & Vasas, 2009; Gouveia et al., 2021). The omnivory index would have posed another interesting index. Like connectivity, it is a measure given for an entire network and rates its stability. Having omnivorous species in a network is important because it makes it possible for the predators to substitute for one another when one species is declining or extinct (Libralato, 2008).

An alternative to using networkx would have been the igraph package. Networkx was chosen as tool for analysis because it provides straight forward functions for different centrality indices. While it is also possible to calculate these indices using igraph, it requires a bit more knowledge of the library functions. However, for certain indices like the omnivory index or the keystone index, there are no functions in networkx either. For more thorough studies, the indices would have to be calculated using mathematical formulas directly either way. An advantage of igraph is that it is compatible with R and C (Igraph Documentation, 2023). For visualising, bokeh has been very useful, however there is no way to depict directed edges using this library. Alternatives would have been using networkx or igraph which would have only been able to produce non interactive plots with very little visual customization options. Pyvis or gravis would have presented interesting approaches, but pyvis has shown errors displaying the graph correctly in jupyter notebook and gravis had limited options. Plotly was the most promising library apart from bokeh, with even more options to choose from. Bokeh was ultimately chosen for being less complex and easier to handle. A problem with using bokeh with networkx was that the locations for the nodes were different every time a graph was generated. Thus, a comparison of the indices is harder to visualise effectively. This happened because no fixed x- and y-coordinates were assigned to the nodes, but the layout was loaded from networkx spring layout.

5 Bibliography

Baeta, A., Niquil, N., Marques, J. C., & Patrício, J. (2011). Modelling the effects of
     eutrophication, mitigation measures and an extreme flood event on estuarine benthic food
     webs. Ecological Modelling, 222(6), 1209–1221.
     https://doi.org/10.1016/j.ecolmodel.2010.12.010

Christensen, V., Walters, C., & Pauly, D. (2005). Ecopath with Ecosim: A User’s Guide.
     Fisheries Centre, University of British Columbia, Vancouver, Canada. www.ecopath.org

Fedor, A., & Vasas, V. (2009). The robustness of keystone indices in food webs. Journal of
     Theoretical Biology, 260(3), 372–378. https://doi.org/10.1016/j.jtbi.2009.07.003

Gouveia, C., Móréh, Á., & Jordán, F. (2021). Combining centrality indices: Maximizing the
     predictability of keystone species in food webs. Ecological Indicators, 126, 107617.
     https://doi.org/10.1016/j.ecolind.2021.107617

Hagberg, A. A., Schult, D. A., & Swart, P. J. (2008). Exploring Network Structure, Dynamics,
     and Function using NetworkX. Proceedings of the 7th Python in Science Conference
     (SciPy2008).

Hui, D. (2012). Food Web: Concept and Applications. Nature Education Knowledge.
     https://www.nature.com/scitable/knowledge/library/food-web-concept-and-applications-
     84077181/

Igraph documentation. (2023). https://igraph.readthedocs.io/en/stable/index.html

Jordán, F., Benedek, Z., & Podani, J. (2007). Quantifying positional importance in food webs:
     A comparison of centrality indices. Ecological Modelling, 205(1), 270–275.
     https://doi.org/10.1016/j.ecolmodel.2007.02.032

Libralato, S. (2008). System Omnivory Index. Encyclopedia of Ecology.
     https://doi.org/10.1016/B978-0-12-409548-9.00605-9

Ortiz, M., Levins, R., Campos, L., Berrios, F., Campos, F., Jordán, F., Hermosillo, B.,
     Gonzalez, J., & Rodriguez, F. (2013). Identifying keystone trophic groups in benthic
     ecosystems: Implications for fisheries management. Ecological Indicators, 25, 133–140.
     https://doi.org/10.1016/j.ecolind.2012.08.020

Pimm, S. L., Lawton, J. H., & Cohen, J. E. (1991). Food web patterns and their consequences.
     Nature, 350(6320), 669–674. https://doi.org/10.1038/350669a0

Rocchi, M., Scotti, M., Micheli, F., & Bodini, A. (2017). Key species and impact of fishery
     through food web analysis: A case study from Baja California Sur, Mexico. Journal of
     Marine Systems, 165, 92–102. https://doi.org/10.1016/j.jmarsys.2016.10.003

Solé, R. V., & Montoya, J. M. (2001). Complexity and fragility in ecological networks.
     Proceedings of the Royal Society B: Biological Sciences, 268(1480), 2039–2045.
     https://doi.org/10.1098/rspb.2001.1767

Vindigni, A. (2022). Web of life tutorial. https://www.web-of-life.es/tutorial/

Wallach, A. D., Dekker, A. H., Lurgi, M., Montoya, J. M., Fordham, D. A., & Ritchie, E. G.
     (2017). Trophic cascades in 3D: Network analysis reveals how apex predators structure
     ecosystems. Methods in Ecology and Evolution, 8(1), 135–142.
     https://doi.org/10.1111/2041-210X.12663