# Exploring Lightweight Federated Learning for Distributed Load Forecasting

Abhishek Duttagupta, Jin Zhao, Shanker Shreejith  Electronic & Electrical
Engineering
Trinity College Dublin, Ireland
Email<EMAIL_ADDRESS>
###### Abstract

Federated Learning (FL) is a distributed learning scheme that enables deep
learning to be applied to sensitive data streams and applications in a
privacy-preserving manner. This paper focuses on the use of FL for analyzing
smart energy meter data with the aim to achieve comparable accuracy to state-
of-the-art methods for load forecasting while ensuring the privacy of
individual meter data. We show that with a lightweight fully connected deep
neural network, we are able to achieve forecasting accuracy comparable to
existing schemes, both at each meter source and at the aggregator, by
utilising the FL framework. The use of lightweight models further reduces the
energy and resource consumption caused by complex deep-learning models, making
this approach ideally suited for deployment across resource-constrained smart
meter systems. With our proposed lightweight model, we are able to achieve an
overall average load forecasting RMSE of 0.17, with the model having a
negligible energy overhead of 50 mWh when performing training and inference on
an Arduino Uno platform.

Index Terms— Federated learning, deep neural networks, Non - i.i.d
distribution, data heterogeneity

##  I Introduction

Electricity demand forecasting is a crucial function in the energy industry to
allow both producers and distribution system operators to optimize the
generation and distribution of electricity in real-time [1]. It aims to
balance the supply and demand curve and maintain the stability of the electric
grid across the entire network of producers and consumers. Smart energy meters
[2] play a vital role in load forecasting by integrating communication
capabilities on top of their ability to monitor energy consumption accurately
[3]. Smart meters transmit demand data from consumers at specified intervals
to distribution system operators, allowing them to observe the data in real
time, make predictions on future demand, and adapt the grid to respond to this
expected demand. This data-driven approach has enabled electric grid operators
to optimize energy generation and storage to adapt to load peaks from the
supply-side, and at the same time achieve demand-side load balancing through
novel schemes like dynamic pricing tariff rates to achieve supply-demand
equilibrium [4].

Researchers have proposed different algorithms for load forecasting that
primarily enables dynamic energy pricing and usage optimization [5]. In [4],
the authors apply gradient-boosted regression trees on the Spanish Energy
dataset to minimize prediction error for a day-ahead forecast, showing the
effectiveness of ML models in load forecasting. However, load forecasting is a
challenging task due to the stochastic nature of customer usage profiles. In
[6], the authors show that the inclusion of spatial and temporal parameters
can have a significant impact on the accuracy of load prediction. Area-based
clustering techniques have also been explored to group customers with similar
electric energy consumption patterns in a spatiotemporal window to improve the
accuracy of load forecasting [7, 8].

Hybrid models combining classical and deep learning methods have also
demonstrated promising results in load forecasting. In [9], the authors
developed a hybrid model utilizing ARIMA, logistic regression, and deep neural
networks to predict peak load days with a 70 % accuracy. Other studies
reported in the literature have proposed models based on Convolutional neural
networks (CNNs) with fuzzy time series, achieving superior results compared to
traditional models [10, 11]. Using complex models such as a standalone Long
Short Term Memory (LSTM) model [12], a combination of LSTM and CNN [13], and
Gated Re-current Units (GRU) combined with CNN [14] have also shown to achieve
similar prediction performance. The use of rolling updates and Bi-LSTMs was
shown to reduce the computational time over other deep learning methods, while
also improving the accuracy of load forecasting [15].

While most deep learning models have been primarily focused on improving
forecasting performance, privacy concerns have risen with the centralized
learning approach used by load forecasting schemes. In the case of centralized
models, data from smart meters of consumers is monitored and accumulated at
locations within the distribution grid, which is used by the model to predict
peak loads or other learned parameters about the demand in specific clusters.
It has been shown that non-intrusive analysis of granular load curves can
provide information about appliances, usage patterns and hence personal
details about consumers [16, 17]. Privacy-preserving machine learning (ML)
through federated learning (FL) is a promising scheme to address the
challenges of centralized learning on sensitive data [17]. With FL, individual
models can be trained collaboratively and locally, removing the need to
accumulate and share sensitive data across the network. While FL has shown to
be effective in other domains [18, 19], its application to load forecasting is
very limited. In [20], the authors propose a combination of LSTM and dense
neural network within an FL framework for load forecasting and utilize
clustering to improve the model’s prediction performance. However, the root
mean squared error (RMSE) measure of the federated approach was higher than
the centralized learning technique. The work in [21] achieved better
prediction performance using a global generalized model aggregated from
locally trained LSTM models at the consumer end, and by incorporating
spatiotemporal factors into the framework. However, the complexity of local
and global models as well as their high energy consumption is a barrier to
their widespread deployment in constrained devices such as smart meters.

In this work, we apply an FL framework for short-term load forecasting using
the openly available smart meter energy consumption dataset for London
households [22]. The key objective is to show that comparable accuracy to
state-of-the-art methods in load prediction can be achieved using lightweight
deep neural network (DNN) models within an FL framework to restrict the need
for granular data sharing and aggregation across the network. To replicate
real-world non-independent and identical distribution (i.i.d.) conditions, we
consider limited random connections for each round of FL, simulating scenarios
where not all devices are participating in the learning process. The main
contributions of this paper are summarized below.

  1. 1.

A lightweight feed-forward model within an FL framework to achieve identical
and/or superior short-term load forecasting performance, even while
considering limited device connections within the cluster (non i.i.d
condition).

  2. 2.

Combining the FL flow with clustering to achieve a better global
representation within each cluster, leading to improved forecasting
performance.

  3. 3.

Evaluate our approach and demonstrate that lightweight models are sufficient
to achieve similar load forecasting performance when compared to other complex
models.

  4. 4.

Quantify the overheads of the proposed lightweight models on a lightweight
microcontroller that mimics an actual smart meter device.

The rest of this paper is organized as follows: Section II discusses the
proposed federated learning framework, and Section III illustrates the
employed methodology in detail; Section IV presents the case study performed
with the results analysis; and conclusions are drawn in Section V.

##  II Light-weighted FL-based Load Forecasting Scheme

###  II-A Federated Learning Algorithm

In a traditional centralized scenario 1, energy consumption data from
participating household smart meters is sent to a cluster head (typically at a
substation level), where the global server performs model training and
inference. In contrast, with an FL setup 2, each household performs its own
edge training of local models based on their individual consumption data.
These updated models are then sent to the global server at the substation for
aggregation, creating the global model. The global model is subsequently
shared with the households for evaluation. The entire cycle keeps repeating,
allowing the local and global model to evolve based on changes in individual
consumers’ consumption and generation.

Figure 1: Traditional centralized setup for load forecasting. Figure 2:
Federated learning setup for load forecasting.

The federated algorithm flow 1 used for our process is derived, from the work
in [23]. The algorithm begins by selecting a set of user devices (N=0,1,2,…,n)
from the total device set (M) and initializing an initial global model (W),
which is then passed to these devices. Training is performed in parallel on
these selected devices using their local data (P), considering batch sizes
(B), learning rate ($\eta$), and local epochs (E). After training, the updated
models (w) from these devices are shared with the substation, where federated
averaging occurs for server aggregation, resulting in a global model (W). The
global model is then shared back with all the devices, which update their
models accordingly, and the process is repeated until convergence is reached.
By Convergence, we refer to the point where no noticeable improvement in the
global loss function (GLF) occurs with increasing federated training rounds,
indicating that the GLF has reached a minimum.

0: Total device set $M$, selected device set $N$, batch size $B$, learning
rate $\eta$, local epochs $E$

0: Global model $W$

1: Initialization: Randomly initialize global model $W$

2: while not converged do

3: for each device $n$ in selected device set $N$ do

4: Retrieve local data $P$ on device $n$

5: Train local model $w_{n}$ on data $P$ for $E$ epochs with batch size $B$
and learning rate $\eta$

6: Share local model $w_{n}$ with central server

7: end for

8: Aggregate models $w_{n}$ from all devices on a central server using
federated averaging

9: Update global model $W$ as the average of all local models $w_{k}$

10: Share global model $W$ with all devices

11: end while

Algorithm 1 Federated learning algorithm

Adopting the representation in [24] for our case, an FL cluster with _n_
clients communicating with the global server(substation for our case), if each
client has a _k_ fixed number of data points in its respective dataset, the
input/output relation for a given local dataset _S_ for a client _i_ can be
represented by the equation eq. 1. Here _x_ and _y_ represent the input and
output vectors, respectively.

$\mathcal{S}_{i}=\\{(x_{i},y_{i})\\}_{i=1}^{k_{i}}$ (1)

The learning process can hence be captured using the model f(.; $\theta$) that
maps an input _x_ to the predicted label _y_. Here, $\theta$ represents the
trainable parameter of the model f. Assuming a loss function of _L_ for this
model capturing the error in the prediction f(x; $\theta$) given the true
label y, the local objective of the client _i_ is given by the equation eq. 2
and the global server objective in an FL setting can be represented by the
equation eq. 3.

$L(\theta;\mathcal{S}_{i})=\frac{1}{k_{i}}\sum_{i=1}^{k_{i}}l(f(x_{i};\theta),y_{i})$
(2)

$\arg\min_{\theta}L(\theta)=\frac{1}{n}\sum_{i=1}^{n}L(\theta;\mathcal{S}_{i})$
(3)

###  II-B Performance Metrics

With errors varying according to household average consumption level, MAE
(calculated using eq. 4) and RMSE (given by eq. 5) are not ideal metrics to
accurately assess short-term load forecasting at small energy consumption
units. Mean absolute percentage error (MAPE, given by eq. 6) provides a more
reliable measure by capturing the absolute percentage deviation between
predicted and actual values. Using MAPE with MAE and RMSE allows for a
comprehensive evaluation and comparison of prediction accuracy independent of
a household’s average consumption.

$MAE=\frac{1}{N}\sum_{i=1}^{N}|X(i)-Y(i)|$ (4)

$RMSE=\sqrt{\frac{1}{N}\sum_{i=1}^{N}[X(i)-Y(i)]^{2}}$ (5)

$MAPE=\frac{1}{N}\sum_{i=1}^{N}\left|\frac{X(i)-Y(i)}{X(i)}\right|\times 100$
(6)

As before, in the case of equations eq. 4 to eq, 6, X represents the actual
value and Y represents the predicted value, with N as the number of data
points.

##  III DISTRIBUTED LOAD FORECASTING METHODOLOGY

For our simulation and testing, we utilized the smart meter energy consumption
dataset from London Households [22]. This dataset comprises energy consumption
data recorded from smart meters installed in 5,567 households across London.
The data covers the period from November 2011 to February 2014, with half-
hourly energy consumption data, containing four columns: a unique tag, tariff
type (standard or dynamic pricing), time-stamp, and half-hourly electric
consumption (in kWh).

First, we examined the heterogeneity of the consumption data and analyzed it
to understand fluctuations in the energy consumption patterns of consumers
across the year and to determine any seasonal variations or long-term trends.
Figure 3 illustrates the monthly total consumption trend for a random
household across the year, indicating that January, February, and December
tend to have higher consumption compared to other months. Similarly, figure 4
shows the daily consumption of one high-consumption consumer and one medium-
consumption consumer, capturing the significant variations and heterogeneity
present in the dataset. We use the average consumption for each household
across the entire dataset to identify outliers with very high or very low
average half-hourly consumption. Based on our observations, households with an
average energy consumption per half-hourly (hh) interval less than 0.09 kWh/hh
(nearly unused) or greater than 1.35 kWh/hh are marked as outliers and
filtered out of the dataset. Eliminating these outliers reduces the number of
households from 5547 to 4672 and removes any bias that may be incurred when
training the model with the outlier data points.

Figure 3: Energy profile yearly trend of a random household Figure 4:
Comparison of daily consumption over a month: high consumption household vs
average consumption household

###  III-A Clustering of smart meter-based energy consumption data

Clustering the training data has shown to provide better generalisability for
deep learning models when used for load prediction [7, 8]. With federated
models, clustering is particularly effective as it allows clients to be
organized into smaller clusters based on similar consumption and their
locality. This also improves the aggregation phase in the substation where
local models from similarly performing households have a better probability of
achieving an unbiased aggregation to arrive at the new global model. For our
approach, we use statistical information about energy consumption to organize
households into 18 groups using K-Means clustering. We use mean energy
consumption, median consumption, total energy consumption, maximum recorded
consumption, and minimum recorded consumption parameters to guide the
clustering process, similar to the approach in [21]. Clustering also takes
into account the area locality for the grouping to ensure that they can be
grouped under a global server that serves this area (e.g., a substation). So,
in our setup, each cluster would indicate an area within London, comprising
households of similar electric consumption levels with a substation (global
server), feeding to all these households. Figure 5 illustrates our constructed
framework up to the cluster level.

Figure 5: Top to Cluster level Simulation Architecture

###  III-B Federated Training and Testing Data sets Creation

The entire dataset, indexed by household tags, is subsequently partitioned
into training and test data. For every household, data from Nov 2011 to Feb
2013 was used for training and Validation purposes, while, the data from Feb
2013 to Feb 2014, was used for testing purposes, thereby achieving a 60:40
split between training and test data.

####  III-B1 Input Window and Output Prediction

To account for the time-series nature of load forecasting, we use a sliding
window scheme to capture the short history of recent consumption data, which
forms the input to the deep learning model. For our experiments, we evaluated
different window sizes using training error at a fixed number of epochs to
arrive at a size of 336, which captures an entire week of consumption at a
half-hour resolution. We infer that by considering a longer time span that
includes weekdays and weekends, we can capture and analyze the trends and
patterns that emerge throughout a complete week.

####  III-B2 Model architecture and TensorFlow Federated Wrappers

For building our training and inference flows (simulation), we utilise the
Python-based TensorFlow-Federated (TFF) library and its associated frameworks.
TFF wrapper converts the windowed input/output data into a Python dictionary
for use with the TFF framework. For our evaluation, we simulate the learning
of each individual household and its aggregation at the cluster using the TFF
libraries. Unique tags associated with each household are used to access
individual data streams and clusters during the training and evaluation phase.

For our model, we chose to evaluate feed-forward neural networks for their
simplicity and computational efficiency for each of the clients. After
evaluating multiple architectures, we settled on a 4-layer model with 2 hidden
layers. The model follows [16,8,4,1] shape with rectified linear unit (ReLU)
activation functions at each layer. We used RMSE and MAE as the loss functions
and metrics for our training and all three (RMSE, MAE and MAPE) for
evaluation. The learning rate was set to 0.01 at each client site. The final
model architecture has 5569 trainable parameters, making it ideal for training
and performing inference on low-cost microcontroller-based computing nodes
(like smart meters) while having sufficient parameters to learn crucial
patterns and trends in the data stream.

##  IV CASE STUDY

###  IV-A Case Settings

We use the following hyper-parameters for training the model: batch size $=$
12, global server learning rate $=$ 1 and number of federated rounds $=$ 20\.
We present a selection of results from individual clusters representing the
spectrum of consumption levels that are available in the dataset in addition
to the aggregate performance across all clusters. For every cluster, two types
of households are selected for evaluating our model performance:

  1. 1.

Households that received all global model updates.

  2. 2.

Households that did not receive some updates or all updates during training.

This differentiation factors in non-i.i.d conditions that may be present in a
real-world condition. Additionally, within the group of households that did
not receive some updates, we further divided them into inter-cluster examples,
including high-consumption and low-consumption examples for observing how the
forecasting accuracy’s varies across different households.

###  IV-B Results

TABLE I: Monthly load forecasting RMSE: comparison between different cluster
08 households

Month | Client Type (RMSE)
---|---
Mod(All) | High(Non-i.i.d) | Low(Non-i.i.d)
January | 0.0636 | 0.2507 | 0.0778
February | 0.0633 | 0.2581 | 0.0665
March | 0.0623 | 0.2875 | 0.0760
April | 0.0544 | 0.1531 | 0.0797
May | 0.0570 | 0.0866 | 0.0740
June | 0.0621 | 0.0583 | 0.0752
July | 0.0638 | 0.1690 | 0.7425
August | 0.0662 | 0.2217 | 0.7119
September | 0.0528 | 0.1417 | 0.0730
October | 0.0611 | 0.1896 | 0.0695
November | 0.0539 | 0.2505 | 0.0684
December | 0.0650 | 0.2484 | 0.0787

Table I presents the RMSE values averaged across different types of households
within a specific cluster from our experiments. It can be observed that the
model offers better (average) prediction performance for households that
received all the updates (i.i.d condition) compared to those that missed out
on some of the updates. Additionally, this demonstrates that in the worse
scenario when a house skip updates, the approach is able to provide reasonably
accurate predictions.

We also examined the performance across festive periods where we observe
significant differences in daily consumption while analyzing the dataset. The
prediction results around a three-day window around Christmas from the 24th of
December to the 26th of December is plotted in figure 6 for different
clusters. The general consumption data is well-aligned with the model’s
prediction, as can be seen from the figure across all clusters. While the
sudden peaks are not accurately predicted by the model, the predicted
consumption is on the lower side than the actual consumption, making it a
reliable baseline for distribution system operators and grid operators.

(b) Cluster 04 (c) Cluster 18 (d) Cluster 07 (e) Cluster 09

(a) Cluster 08

Figure 6: Short-term daily forecasting (24th – 26th of December) across a
selection of clusters

####  IV-B1 Performance compared with Centralized Setup

We evaluate the use of our lightweight model in a purely centralized setup to
see the effectiveness of distributed learning in comparison to a global
learning approach using a model of similar complexity. For the centralized
setup, we assume that all consumption data is made available to the central
server and that all clients participate at a consistent rate (set to 3$\times$
higher) compared to our federated learning scheme. The clustering approach is
maintained in the centralized setup to ensure fairness in both cases. Both
setups use 20 rounds of training - FL using training locally at the clients,
while the centralized model trains for 20 epochs at the central server. Figure
7 shows the average prediction error (RMSE) for the entire year for a single
cluster (cluster 04) across different months. The results show that the
federated model achieves better performance at the same model complexity.

Figure 7: RMSE of centralized vs federated prediction for cluster 04

####  IV-B2 Performance comparison against other research works

We compare our load forecasting model with the federated architecture used in
[21], using average monthly RMSE across various clusters. We include
households with high, low and moderate electric energy consumption. For our
simulation, the number of participating clients was set to be less than 50% of
the clients used in the competing method [21] for each cluster. Although an
increase in the number of participating nodes does not guarantee improvement
in prediction performance, it should potentially lead to a better global model
based on the consumption profile of the additionally included nodes.
Additionally, the work in [21] uses an LSTM-based model with 7505 parameters
compared to our 4-layer feed-forward model with only 5569 learnable parameters
(and thus lower communication costs during the model aggregation phase).

TABLE II: Comparison of average RMSE of all clusters between our 4-layer feed-
forward model and the LSTM-based FL model in [21]

Month | Our Proposed Model (RMSE) | M. Savi Model (RMSE)
---|---|---
January | 0.2053 | 0.1463
February | 0.1626 | 0.153
March | 0.2176 | 0.1535
April | 0.1800 | 0.1345
May | 0.1657 | 0.1304
June | 0.1434 | 0.1259
July | 0.1685 | 0.1057
August | 0.1526 | 0.1106
September | 0.1581 | 0.1302
October | 0.1554 | 0.1335
November | 0.1729 | 0.1471
December | 0.1728 | 0.1328

The table II) shows the results of our simulation capturing the average
monthly RMSE of the models across all clusters over the entire year in
comparison with the results reported in [21]. Despite using 50% fewer
participating nodes in each cluster and 20% fewer computational requirements
for the local model, our model is able to achieve an average RMSE of 0.17
compared to 0.14 of the LSTM model. While RMSE does not fully capture the
prediction performance in this case, we believe that the trade-off is in line
with the lower computational complexity, communication overhead and
participation ratio in the network.

We also compare our results against the work in [25], where the model uses a
participation ratio of 36% compared to the 10% in our method. Similar to the
previous work, the work in [25] also uses an LSTM-based model with customized
local model improvement at each round, with only 5 federated rounds compared
to our 20. Our results, shown in table III, show that our approach achieves
significantly better average prediction accuracy (in terms of MAPE) compared
to the results reported in [25]. Additionally, we observe that reducing the
number of federated rounds from 20 to 5 increases the average error across all
clusters to increase by nearly 4% in our case, which is still better than the
results achieved by [25].

TABLE III: Comparison of Participation Ratio, Number of Rounds, and Average
MAPE

Model | Household Ratio | Fed. Rounds | Avg MAPE
---|---|---|---
Taik’s Model | 36% | 5 | 34.14%
Our Proposed Model | 10.50% | 20 | 22.01%

####  IV-B3 Energy consumption of the model on an IoT platform

To quantify the energy overhead in performing model training and inference at
each client (smart meter), we implement the training and inference on an
Arduino Uno R4 WiFi platform, that uses an 8-bit processor. The model was
coded using Ardiuno C and compiled to the device using the standard compiler
flow. For evaluation, we simulate the condition that the model is training on
the data for one day and predicted the half-hourly consumption for the
following day, with the output transferred via serial port to a standard
laptop. To measure the energy consumed, we monitor the power consumed by the
Arduino device using a calibrated USB energy monitor and record the
consumption of the model under idle, training and prediction phases. We
observed that the model consumes 50 mWh (milliwatt-hour) when averaged over
the entire day on top of the idle power consumption of 440 mWh, when
performing training and inference, resulting in an overhead of 11% on the
Arduino platform. The result indicates that integration of our model will not
incur significant energy consumption overhead in a similarly capable (or
higher) smart meter hardware platform, making our approach ideally suited for
this application.

##  V Conclusion

In this work, we propose a federated learning-based approach for distributed
load forecasting that uses a lightweight fully connected neural network at its
core and test it on a real-world openly available energy consumption dataset.
By combining this lightweight model with clustering and randomized client
selection, we emulate real-world conditions in an electric distribution grid.
Our approach learns locally at each client’s smart meter with only the model
updates passed up to the central entity, reducing the privacy risks associated
with granular aggregation of energy consumption data. Our simulations show
that the light-weight model when combined with clustering is able to generate
accurate load forecasting, and provides comparable results against competing
FL-based load prediction approaches, trading off inference accuracy
($\approx$21.4% higher RMSE) with model complexity ($\approx$25.8% lower
parameters) in one case, while achieving $\approx$35.6% better aggregate MAPE
in another. We also show that training and inference using our model on an
Arduino platform incurs an energy overhead of only 50 mWh when averaged over
multiple cycles, which is a fraction of the standalone power consumption for
the platform, making it ideally suited for load forecasting in a decentralized
setting. In the future, we aim to show that decentralized learning can aid in
improving the smart-grids stability and resilence without requiring
aggregation of privacy-sensitive data.

## References

  * [1] R. Weron, “Electricity price forecasting: A review of the state-of-the-art with a look into the future,” _International journal of forecasting_ , vol. 30, no. 4, pp. 1030–1081, 2014.
  * [2] A. Arif, M. Al-Hussain, N. Al-Mutairi, E. Al-Ammar, Y. Khan, and N. Malik, “Experimental study and design of smart energy meter for the smart grid,” in _2013 International Renewable and Sustainable Energy Conference (IRSEC)_. IEEE, 2013, pp. 515–520.
  * [3] E. Spanò, L. Niccolini, S. Di Pascoli, and G. Iannacconeluca, “Last-meter smart grid embedded in an internet-of-things platform,” _IEEE Transactions on smart grid_ , vol. 6, no. 1, pp. 468–476, 2014.
  * [4] G. Díaz, J. Coto, and J. Gómez-Aleixandre, “Prediction and explanation of the formation of the spanish day-ahead electricity price through machine learning regression,” _Applied Energy_ , vol. 239, pp. 610–625, 2019.
  * [5] A. R. Khan, A. Mahmood, A. Safdar, Z. A. Khan, and N. A. Khan, “Load forecasting, dynamic pricing and dsm in smart grid: A review,” _Renewable and Sustainable Energy Reviews_ , vol. 54, pp. 1311–1322, 2016\.
  * [6] A. Tascikaraoglu and B. M. Sanandaji, “Short-term residential electric load forecasting: A compressive spatio-temporal approach,” _Energy and Buildings_ , vol. 111, pp. 380–392, 2016.
  * [7] K. Bandara, C. Bergmeir, and S. Smyl, “Forecasting across time series databases using recurrent neural networks on groups of similar series: A clustering approach,” _Expert systems with applications_ , vol. 140, p. 112896, 2020.
  * [8] F. L. Quilumba, W.-J. Lee, H. Huang, D. Y. Wang, and R. L. Szabados, “Using smart meter data to improve the accuracy of intraday load forecasting considering customer behavior similarities,” _IEEE Transactions on Smart Grid_ , vol. 6, no. 2, pp. 911–918, 2014.
  * [9] H. Saxena, O. Aponte, and K. T. McConky, “A hybrid machine learning model for forecasting a billing period’s peak electric load days,” _International Journal of Forecasting_ , vol. 35, no. 4, pp. 1288–1303, 2019\.
  * [10] G. Sideratos, A. Ikonomopoulos, and N. D. Hatziargyriou, “A novel fuzzy-based ensemble model for load forecasting using hybrid deep neural networks,” _Electric Power Systems Research_ , vol. 178, p. 106025, 2020.
  * [11] H. J. Sadaei, P. C. d. L. e Silva, F. G. Guimarães, and M. H. Lee, “Short-term load forecasting by using a combined method of convolutional neural networks and fuzzy time series,” _Energy_ , vol. 175, pp. 365–377, 2019.
  * [12] G. Memarzadeh and F. Keynia, “Short-term electricity load and price forecasting by a new optimal lstm-nn based prediction algorithm,” _Electric Power Systems Research_ , vol. 192, p. 106995, 2021.
  * [13] X. Guo, Q. Zhao, D. Zheng, Y. Ning, and Y. Gao, “A short-term load forecasting model of multi-scale cnn-lstm hybrid neural network considering the real-time electricity price,” _Energy Reports_ , vol. 6, pp. 1046–1053, 2020.
  * [14] M. Sajjad, Z. A. Khan, A. Ullah, T. Hussain, W. Ullah, M. Y. Lee, and S. W. Baik, “A novel cnn-gru-based hybrid approach for short-term residential load forecasting,” _Ieee Access_ , vol. 8, pp. 143 759–143 768, 2020.
  * [15] S. Wang, X. Wang, S. Wang, and D. Wang, “Bi-directional long short-term memory method based on attention mechanism and rolling update for short-term load forecasting,” _International Journal of Electrical Power & Energy Systems_, vol. 109, pp. 470–479, 2019.
  * [16] G. Eibl and D. Engel, “Influence of data granularity on smart meter privacy,” _IEEE Transactions on Smart Grid_ , vol. 6, no. 2, pp. 930–939, 2015.
  * [17] L. AlAbdulkarim and Z. Lukszo, “Impact of privacy concerns on consumers’ acceptance of smart metering in the netherlands,” in _2011 International Conference on Networking, Sensing and Control_ , 2011, pp. 287–292.
  * [18] B. Hu, Y. Gao, L. Liu, and H. Ma, “Federated region-learning: An edge computing based framework for urban environment sensing,” in _2018 IEEE Global Communications Conference (GLOBECOM)_. IEEE, 2018, pp. 1–7.
  * [19] S. R. Pfohl, A. M. Dai, and K. Heller, “Federated and differentially private learning for electronic health records,” _arXiv preprint arXiv:1911.05861_ , 2019.
  * [20] N. Gholizadeh and P. Musilek, “Federated learning with hyperparameter-based clustering for electrical load forecasting,” _Internet of Things_ , vol. 17, p. 100470, 2022. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S2542660521001104
  * [21] M. Savi and F. Olivadese, “Short-term energy consumption forecasting at the edge: A federated learning approach,” _IEEE Access_ , vol. 9, pp. 95 949–95 969, 2021.
  * [22] UK Power Networks: SmartMeter Energy Consumption Data in London Households, “https://data.london.gov.uk/dataset/smartmeter-energy-use-data-in-london-households,” 2014\.
  * [23] B. McMahan, E. Moore, D. Ramage, S. Hampson, and B. A. y Arcas, “Communication-efficient learning of deep networks from decentralized data,” in _Artificial intelligence and statistics_. PMLR, 2017, pp. 1273–1282.
  * [24] Z. Shen, J. Ye, A. Kang, H. Hassani, and R. Shokri, “Share your representation only: Guaranteed improvement of the privacy-utility tradeoff in federated learning,” in _The Eleventh International Conference on Learning Representations_ , 2023.
  * [25] A. Taïk and S. Cherkaoui, “Electrical load forecasting using edge computing and federated learning,” in _ICC 2020-2020 IEEE international conference on communications (ICC)_. IEEE, 2020, pp. 1–6.


# Hybrid structural arrangements mediate stability and feasibility in
mutualistic networks

Aniello Lampo  Universitat Rovira i Virgili, Tarragona, Catalonia, Spain
María J. Palazzi  Internet Interdisciplinary Institute (IN3), Universitat
Oberta de Catalunya, Barcelona, Catalonia, Spain  Javier Borge-Holthoefer
Internet Interdisciplinary Institute (IN3), Universitat Oberta de Catalunya,
Barcelona, Catalonia, Spain  Albert Solé-Ribalta  Internet Interdisciplinary
Institute (IN3), Universitat Oberta de Catalunya, Barcelona, Catalonia, Spain
URPP Social Networks, University of Zurich, Zurich, Switzerland

###### Abstract

The debate about ecosystems persistence, namely how ecological communities
sustain a large number of species, holds considerable importance since the
middle of the last century. It has been recognized that the persistence of a
system requires stability, but also feasibility, defined as the capability to
sustain positive abundances for all biotas. So far, the relationship between
the emergence of ecological variables and the species interaction networks has
been studied by restricting to single structural patterns, such as nestedness
or modularity. Focusing on mutualistic communities, the current paper
addresses the role of the interaction architecture in the emergence and
maintenance of both properties. Specifically, we propose a new approach
relying on hybrid architectural configurations, motivated by the intuition
that the co-occurrence of multiple variables cannot be properly grasped by
considering different network configurations distinctly. We launch in-block
nestedness in the field of ecology, compound by disjoint subsets of species
(modules) with internal nested character, and prove that it trades off
stability and feasibility. Remarkably, we analyze a large amount of empirical
communities and find that a relevant fraction of them exhibits a marked in-
block nested structure. We elaborate on the implications of these results,
arguing that they provide useful insights about the key properties ruling
community assembly.

##  I INTRODUCTION

Ecological systems are constituted by a high number of species that share a
habitat, and interact with each other in varied ways –predation, mutualism,
competition, etc. Interactions among species play a crucial role in the
dynamics and organization of communities. They shape the ongoing evolutionary
process and represent a key component of several properties, for example, the
maintenance of diversity, the coexistence and abundance of species, as well as
the stability of the whole ecosystem. Herein, network theory has contributed
important advances in the understanding of how species interactions affect
these variables.

Previous studies pointed out that interaction among biotas are not randomly
assembled, but manifest clear architectural patterns. For instance, Bascompte
et al. Bascompte  _et al._ (2003) analyzed 52 plant-animals mutualistic
communities finding that they are organized in nested arrangements Atmar and
Patterson (1993); Mariani  _et al._ (2019), i.e. the specialist species
interact only with subsets of those interacting with the more generalists.
Similarly, Olesen et al. Olesen  _et al._ (2007) showed as well that
pollination networks can exhibit a modular character Newman (2006), with
weakly interlinked groups of species (modules), which exhibit strong internal
connectivity. Besides these two common examples of interaction configuration,
there exist several other types that remain largely understudied Lewinsohn
_et al._ (2006).

The evidence that real ecosystems networks can be strongly structured
triggered a long series of research investigating the relationship between
these structural arrangements and several dynamical property of the ecosystem.
Much effort was devoted to explore the effect of interaction architecture
Okuyama and Holland (2008); Thébault and Fontaine (2010); Allesina and Tang
(2012); Stouffer and Bascompte (2011); Staniczenko  _et al._ (2013); Allesina
_et al._ (2015); Grilli  _et al._ (2016) on stability. In particular, in
Allesina and Tang (2012); Staniczenko  _et al._ (2013) the authors pointed out
that nested networks are less likely to be stable, while the several other
works discussed the stabilizing effects associated to modularity Stouffer and
Bascompte (2011); Allesina  _et al._ (2015); Grilli  _et al._ (2016). Leaving
stability aside, Bastolla et al. Bastolla  _et al._ (2009) proved that
nestedness promotes diversity. Further, nested arrangements have been shown to
positively correlate with ecosystems feasibility Saavedra  _et al._ (2015), a
quantity related to species coexistence, or be a result of optimizing species
abundances Suweis  _et al._ (2013).

Taking a general overlook at the state-of-the-art, it is easy to realize that,
so far, the structure-dynamics interplay has been centered on one-to-one
mappings, examining the correspondence of a given dynamical property to a
single architectural pattern, or viceversa. This necessarily leads to a
description of the resulting behavior where a particular property predominates
over the remaining ones. For instance, the analysis in terms of modularity
portrays a picture of a community in which stability is promoted over
diversity Staniczenko  _et al._ (2013) or, similarly, nestedness benefits
feasibility Saavedra  _et al._ (2015). However, this reductionistic view may
not fit with the more complex dynamics of natural systems, which seem to
evolve by jointly optimizing several ecological variables Tilman and Downing
(1994); McCann (2000). Note that several structural arrangements, offering
complementary beneficial aspects to an ecosystem, may not be structurally
compatible with each other Palazzi  _et al._ (2019a), pushing this path to a
dead-end. Limited to a single dynamical property, multilayer or multiplex
networks Gracia-Lázaro  _et al._ (2018); Stella  _et al._ (2018); Pilosof  _et
al._ (2017); Kéfi  _et al._ (2016); García-Callejas  _et al._ (2018);
Hutchinson  _et al._ (2019) may address some aspects of the problem. Here, we
seek a different research direction, exploring hybrid structural arrangements.

The study of how nested and modular network arrangements may entangle has
become, during the last years, a central issue empirically Olesen  _et al._
(2007); Fortuna  _et al._ (2010); Flores  _et al._ (2011, 2013); Beckett and
Williams (2013) and theoretically Solé-Ribalta  _et al._ (2018); Palazzi  _et
al._ (2019b); Cai  _et al._ (2020); Pinheiro  _et al._ (2019). Particular
attention has been payed to the coexistence of both patterns because, as
discussed, they seem arise as a consequence of different evolutive mechanisms.
Here, in-block nestedness (IBN) stands out as a potential hybrid network
configuration, consisting of blocks whose internal connectivity is arranged in
a nested manner. Such compound architecture was first described in Lewinsohn
_et al._ (2006), and structurally embodies properties of both modularity and
nestedness Solé-Ribalta  _et al._ (2018); Palazzi  _et al._ (2019b).

The main aim of our work is to show that hybrid structures are crucial to
approach the ecosystem functioning. These structures may offer an interesting
balance between different dynamical properties incompatible a priori, as
implied from the current state-of-the-art. We will consider the case of IBN,
showing that it trades-off stability and feasibility and constitutes an
optimal pattern to ensure the persistence of ecosystems. Such result is
obtained for an ensemble of synthetic networks of different size and
connectance, and holds both in a pure mutualistic context and in presence of
intra-species competition. Furthermore, we analyze almost one hundred and
fifty mutualistic empirical networks, including anemone-fish, plant-pollinator
and seed-dispersal communities, and find that a high fraction of them exhibits
an appreciable degree of IBN –often higher than any other structural
arrangement. This has important consequences on the understanding of the
organization of mutualistic communities, suggesting for many systems the
emergence of hybrid interaction patterns in the attempt to preserve both
stability and diversity, rather than promoting one over the other.

Figure 1: Examples of adjacency matrixes associated to nested (left) modular
(middle) and IBN (right) networks presented in Solé-Ribalta  _et al._ (2018).
For nested networks, the adjacency matrix manifests a triangular structure,
while in the modular case it is divided in blocks, exhibiting a high internal
connectivity and a low external one. IBN matrixes are divided in blocks with
an internal nested structure. Vertical and horizontal lines are a visual guide
to show the existing blocks.

##  II METHODS

###  II.1 Generalized Lotka-Volterra dynamcs

An ecological community may be modeled as a network in which nodes represent
the species, while edges encode the interaction type and strength among them:
a link is turned on between nodes $i$ and $j$ if the related species interact.
On the basis of this interaction network, species abundances can be described
by a set of time-dependent functions $x_{i}(t)$, and their temporal evolution
is commonly studied by recalling the generalized Lotka-Volterra equation

$\frac{dx_{i}}{dt}=x_{i}\left(r_{i}-\sum_{j}{M}_{ij}x_{j}\right),\quad i,j\leq
S,$ (1)

where $S$ corresponds to the species number. The parameters $r_{i}$ represent
the intrinsic self-growth coefficients, ruling the dynamics of the $i$-species
when interspecies interaction is dropped out. The information concerning this
interaction is embodied in matrix $M$, properly endowed by weights, which
corresponds to the adjacency matrix of the interaction network.

In the present work we focus on bipartite mutualistic networks, in which nodes
are divided in two groups $A$ and $P$. Accordingly, the interaction between
different sets species is assumed to be mutualistic, while within species
corresponding to the same set only competition occurs. This captures the
evolution of several bipartite communities of plants and animals, such as
plant-pollinators or seed-dispersal, where species belonging to different sets
cooperate, while compete with those of the same group to get resources. In
this case, the adjacency matrix of the interaction has a particular shape,

$M=\begin{pmatrix}\Omega_{AA}&-\Gamma_{AP}\\\
-\Gamma_{PA}&\Omega_{PP}\end{pmatrix},$ (2)

where the block $\Omega^{AA}$ $\left(\Omega^{PP}\right)$ stands as the
competitive interactions among species belonging to the set $A$ ($P$), and the
block $\Gamma^{AP}$ describes the mutualistic interaction between species
corresponding to different sets. For the sake of simplicity, during the rest
of the manuscript we will assume that both the set $A$ and $P$ contains the
same number of species $S$.

The main aim of the work is to understand how the joint behavior of stability
and feasibility is affected by the architecture of mutualistic interactions.
Thus, the off-diagonal block $\Gamma^{AP}$ (exhibiting the same structure of
$\Gamma^{PA}$) is the relevant structure for our analysis. This is generated
by means of the model introduced in Solé-Ribalta  _et al._ (2018) which
permits to span nested, modular, and IBN structures in terms of four
parameters: the number of communities $B\in\left[1,\infty\right)$, the noise
regarding the existence of interactions outside the communities
$\mu\in\left[0,1\right]$, how nested the interaction within communities is
$p\in\left[0,1\right]$, and a shape parameter controlling the slightness of
the nested structure $\xi\in\left[1,\infty\right]$ (see Fig. 2 in Palazzi  _et
al._ (2019b)). More details regarding the synthetic network generation model
are provided in Appendix A.

The diagonal blocks, following the state-of-the-art Saavedra  _et al._ (2015);
Suweis  _et al._ (2013), are treated as unstructured interactions. The weights
of both the competitive and mutualistic interaction are treated in the mean
field interaction:

$\displaystyle\Gamma^{AP}_{ij}=\langle\Gamma^{AP}_{ij}\rangle+\sigma^{\gamma}_{ij}\equiv\gamma+\sigma^{\gamma}_{ij},\quad\sigma^{\gamma}_{ij}\ll\gamma,$
(3)

and similar for the other blocks. In particular, we assume
$\langle\Omega^{AA}_{ij}\rangle=\langle\Omega^{PP}_{ij}\rangle\equiv\omega$,
namely the average value of the intraspecies competition is the same for the
two sets. In agreement with the previous literature Allesina and Tang (2012);
Saavedra  _et al._ (2015) their diagonal is set to one, i.e. $M_{ii}=1$: these
elements correspond to the species auto-interaction, namely they specify the
behavior of a species due to internal interaction among different specimens.

###  II.2 Stability

Stability is evaluated by looking into the Jacobian matrix of the generalized
Lotka-Volterra model, Eq. 1,

$J_{ij}\equiv\left(\frac{\partial\dot{x}_{i}}{\partial
x_{j}}\right)_{\textbf{x}=\textbf{x}^{*}},$ (4)

where $\textbf{x}^{*}$ represents the stationary state of Eq. (1), defined by
the condition $\dot{x}^{*}_{i}=0,$ leading to

$\textbf{x}^{*}=M^{-1}\textbf{r},$ (5)

and describing the species abundances at equilibrium. Replacing Eq. (1) into
Eq. (4) and recalling Eq. (5) it turns

$J_{ij}=-x^{*}_{i}M_{ij}.$ (6)

The system is stable if the maximum real part of the eigenvalues of the
Jacobian matrix is negative, otherwise is said to be unstable. Hence, the
quantity

$\lambda^{*}=-\text{Max}\left[Re\\{\lambda_{J}\\}\right],\quad\lambda_{J}\in\text{Sp}(J)$
(7)

naturally describes the stability of the system.

The Jacobian matrix provides information about the long-time response of the
first-order dynamical behavior of the system, after infinitesimal
perturbations of the stationary abundances. These in general affect the
expression of the Jacobian, as results from Eq. (6), and have to be taken into
account in the eigenvalues calculation. Following the approach presented in
Gibbs  _et al._ (2018) we assume that the abundances are all positive, i.e.
$x^{*}_{i}>0$ $\forall i\leq S$. This is equivalent to sample a proper r
vector in the feasibility domain, namely to tune the values of the self-growth
rates in order to obtain positive abundances for a given $M$ matrix.

###  II.3 Feasibility

The feasibility regards the conditions on the external parameters yielding a
stationary state in which all species abundances are positive, i.e. $x_{i}>0$,
$\forall i\leq 2S$. In the case of the Generalized Lotka-Volerra model, in Eq.
(1), is easy to infer that, once $M$ is fixed, feasibility only depends on the
self-growth vector r. The feasibility space is spanned by the values of this
quantity associated to positive abundances. This, as pioneered by Svirezhev
and Logofet (1983) and then also shown by Saavedra  _et al._ (2015); Grilli
_et al._ (2017), takes the form of a cone, which borders are given by the
column vectors of $M$.

The degree of feasibility is provided by the openness of the cone, i.e. the
related solid angle $\Omega$, that is completely determined by the $M$ matrix.
Apart of its geometric meaning, $\Omega$ may be interpreted, once it has been
properly normalized, as the probability of sampling randomly a r vector
related to positive abundances. This quantity has been calculated in Ribando
(2006) proving that it takes the following form

$\Omega=\frac{1}{(2\pi)^{S/2}\sqrt{\text{det}(\Sigma)}}\int...\int_{\mathbb{R}^{S}\geq
0}e^{-\frac{1}{2}\textbf{x}^{t}\Sigma^{-1}\textbf{x}}d\textbf{x},$ (8)

which is the cumulative function of a multivariate normal distribution with
mean value equal to zero and variance matrix $\Sigma^{-1}=2M^{t}M$.

Following the approach presented in Saavedra  _et al._ (2015), the $\Omega$
function may be calculated in R language by means of the function pmvnorm,
which involves the computation of multivariate normal probabilities with
arbitrary correlation matrices. This is efficiently performed through quasi-
Monte Carlo method, recalling the algorithm presented in Genz and Bretz
(2009). Hereafter, we represent the solid angle on a logarithmic scale
$\mathcal{F}\equiv\log_{10}\left(\Omega\right)$: this is the quantity we will
actually mean by the term feasibility. It turns that $-\infty<\mathcal{F}\leq
0$ because $0\leq\Omega\leq 1$, since it is a normalized distribution. The
higher the value of $\mathcal{F}$ the larger the likelihood of finding the
system in a feasible state.

Figure 2:  Stability and Feasibility values for an ensemble of synthetic
networks with in-block, nested, and modular architecture. networks refer to
mutualistic communities with species number $S=30$ and connectance $C=0.1$.
Plot on the left describes a pure mutualistic context while that on right
includes intraguild competition, with the same parameters values employed in
Saavedra  _et al._ (2015). Each point represents a network: its horizontal
coordinate is given by the feasibility $\mathcal{F}$ introduced in Sec. II.3,
while the vertical one is provided by the quantity $\lambda^{*}$ introduced in
Eq. (7). It is possible to recognize three different clusters corresponding to
different structures: nested and modularity respectively promote feasibility
and stability, while IBN trades-off the two properties. Red daggers are
located at the average values of stability and feasibility, introduced in Eq.
(9) of the structured networks cluster. Colorbar refers to the $\xi$ values of
the nested (left) and in-block (right) networks.

##  III RESULTS

###  III.1 Synthetic networks

The results of the stability and feasibility calculation for nested, in-block,
and modular synthetic networks are presented in Fig. 2. It is possible to see
that modularity is the more likely structure to be stable, while nestedness
constitutes the less one. In this sense, stability defines a hierarchy in
which modularity is on the top and nestedness in the bottom. Feasibility also
discriminates the three architectures and yields a hierarchy where, at odds of
stability, nestedness is on the top and modularity in the bottom. Importantly,
in both situations IBN is in the middle, standing out as a trade-off between
stability and feasibility. Summarizing, we have

$\langle\mathcal{F}_{\mathcal{N}}\rangle>\langle\mathcal{F}_{\mathcal{I}}\rangle>\langle\mathcal{F}_{\mathcal{Q}}\rangle,\quad\quad\langle\lambda^{*}_{\mathcal{N}}\rangle<\langle\lambda^{*}_{\mathcal{I}}\rangle<\langle\lambda^{*}_{\mathcal{Q}}\rangle,$
(9)

where $\langle\mathcal{F}_{\mathcal{X}}\rangle$ indicates the average value of
feasibility over an ensemble of synthetic networks with $\mathcal{X}$
structure. Similarly for $\langle\lambda^{*}_{\mathcal{X}}\rangle$.

This finding agrees with the previous literature, where nestedness was
recognized to hinder stability Allesina and Tang (2012); Staniczenko  _et al._
(2013). From the mathematical point of view, this fact may be understood in
terms of the Gershgoring theorem, linking the maximum real part of the
Jacobian eigenvalues with the matrix average row sum (related to $\xi$
parameter), which is maximized in nested networks because of the presence of
generalist species. Along this line, we argue that blocks structure minimizes
the row sum, while in-block generalists species permits to reach an
intermediate score between nestedness and modularity. This is certainly true
for the pure mutualistic systems but, worth highlighting, our results are
robust for an average competition strength approximately equal to the $10\%$
of the mutualistic one.

Figure 3:  Stability and Feasibility dependence on the connectance for nested,
in-block and modular networks, at different values of the species number. Each
point refers to a set of matrixes of fixed connectance, and is located at the
resulting average value of feasibility (up) or stability (down). Shadow areas
provides an information about the variance. All the plots refer to the
$\gamma=0.13$ and $\omega=0.01$.

The result presented in Fig. 2 refers to a fixed value of connectance and
species number. In Fig. 3 we extend this analysis over a wide range of such
parameters. It is possible to point out that, for the given values of
competition and mutualism, stability hierarchy holds over the whole selected
domain. Surprisingly, the feasibility-architecture relationship is not uniform
and tends to reverse as $C$ and $S$ increase. A general representation of the
validity region of the feasibility hierarchy is portrayed in Fig. 4. We infer
the the trade-off role of IBN is not universal and concerns a finite, albeit
significative, portion of the parameters space, mainly related to small values
of connectance or species number.

Remarkably, as connectance and species number grow, we approach a situation
where modular networks promote feasibility, as well as stability, and stands
out as optimal interaction architecture to preserve ecological persistence. We
have:

$\langle\mathcal{F}_{\mathcal{N}}\rangle>\langle\mathcal{F}_{\mathcal{I}}\rangle>\langle\mathcal{F}_{\mathcal{Q}}\rangle,\quad\quad\langle\lambda^{*}_{\mathcal{N}}\rangle>\langle\lambda^{*}_{\mathcal{I}}\rangle>\langle\lambda^{*}_{\mathcal{Q}}\rangle,$
(10)

The non-uniform behavior of feasibility also concerns its dependence on
connectance. While for small $S,C$ values feasibility increases as a function
of this quantity at the odds of stability, they approximately match the same
behavior in the opposite limit.

Figure 4: Validity of the trade-off property of IBN in the size-connectance
domain, with $\gamma=0.13$ and $\omega=0.01$. Grey area refers to the region
in which the hierarchy presented in Eq. 9 is fulfilled and IBN constitutes a
trade-off between stability and feasibility, where blue one describes the
parameters space where such property does not hold anymore. Scatter plot
depicts the distribution of empirical pollinators and seed-dispersal
mutualistic communities of wol , with color portraying the in-block nestedness
degree in Eq. (B).

###  III.2 Empirical Networks

An important part of the current work is constituted by the comparison between
our theoretical findings and the properties of real mutualistic ecosystems. We
look into pollinators and seed-dispersal communities of the Web of Life
database wol . This contains 186 networks, ranging over a large domain of size
and connectance. In this work, we restrict to parameters space bounded by
$S_{max}=110$ and $C_{max}=0.3$, where are concentrated 129 communities.

In the previous subsection we showed that the relationship between feasibility
and network architecture is not uniform and depends on the size and
connectance. Fig. 4 presents a comparison between the empirical communities
distribution in the $S$-$C$ space, and the validity regime of the hierarchies
in Eq. (9) for synthetic networks. It is possible to see that the majority of
real networks, i.e. the $90.7\%$, falls into the region where IBN trades off
stability and feasibility, representing a profitable interaction architecture
to preserve ecological persistence.

We are interested in the structural properties of the empirical communities we
considered, namely if they are predominantly nested, modular or in-block. A
suitable manner to present this information is given by a ternary plot (or
simplex), see Fig. 5. Here, each ecological community is located with the
three coordinates $f_{\mathcal{N}}$, $f_{\mathcal{Q}}$, and $f_{\mathcal{I}}$,
defined as
$f_{\mathcal{N}}=\mathcal{N}/(\mathcal{N}+\mathcal{Q}+\mathcal{I})$, and
similarly for $f_{\mathcal{Q}}$, $f_{\mathcal{I}}$, with $\mathcal{N}$,
$\mathcal{Q}$, $\mathcal{I}$ providing the degree of nestedness, modularity
and in-block nestedness (see Sec. B). Note that the, by definition, the
quantities $f_{\mathcal{N}}$, $f_{\mathcal{Q}}$, and $f_{\mathcal{I}}$ take
values in the $[0,1)$ range. The simplex may be divided into dominance
regions, bounded by the three angle bisectors, associated to the several
architectures. Note the part of the ternary are necessarily empty, because of
the constraints mentioned in Appendix B and carefully discussed in Palazzi
_et al._ (2019b).

Fig. 5 unveils that most of the ecological communities fall into the IBN area.
Here, networks result to be relatively small, while those with a notable size
tend to be located in the modularity region. Interestingly, as discussed for
Fig. 4, small communities fulfill the hierarchies in Eq. (9) and IBN exhibits
the trade-off property. Our empirical analysis suggests so that mutualistic
ecosystems develop a IBN architecture to maintain an equilibrium between
stability and feasibility, in order to ensure persistence.

Figure 5:  Distribution of pollinator and seed-dispersal empirical networks
selected by wol , across a ternary plot. The colourbar indicates the in-block
nestedness degree of the networks, evaluated according the measure introduced
in Solé-Ribalta  _et al._ (2018). The point size provides an information of
the species number of the corresponding community.

##  IV DISCUSSION

Stability and feasibility have been mapped, in the past, onto single
preferential architectural patterns. However, the interest in studying
abundant, diverse and stable ecological communities calls for a new vision of
species interaction patterns, able to reproduce the balance of several
properties, rather than the predominance of a single one. With this objective,
we leave the traditional approach focused on a single architectural pattern,
and look into compound structures. We consider in-block nested configuration,
where nestedness and modularity interfere at a network’s mesoscale, and show
that it provides an optimal middle-term between stability and feasibility.
Moreover, the analysis of real mutualistic networks proves that a large
fraction of them actually exhibits in-block nested assemblage.

This result has important consequences on the general understanding of the
mechanisms ruling the organization of mutualistic communities. So far, the
question has been addressed as finding the key property shaping ecosystems
assemblage. For instance, in Saavedra  _et al._ (2015) authors conclude that
communities promote feasibility over stability, but not both at the same time.
This follows from an approach to the study of interaction involving only one
possible structure (nestedness). The introduction of an intertwined
perspective for interaction architecture permits to revisit and deepen this
finding. Particularly, the emergence of in-block nestedness in real
communities, associated to a trade-off between stability and feasibility,
implies that the fundamental criterium underlying the assembly process is the
equilibrium between those properties, rather than the predominance of one of
them.

The mediating role of in-block nestedness depends in general on the
connectance and species number. As these quantities become larger, we detect
an unexpected effect consisting in a reversing of the feasibility hierarchy.
In this regime both stability and feasibility result to be promoted by
modularity, which stands out as the optimal pattern for ecosystems assemblage.
The transition from an in-block to modular optimal pattern may be interpreted
as a witness of the non-universal behavior of ecosystems, which cannot be
treated regardless of its structural properties, such as size. Particularly,
this result paves the way to future studies aimed to look for the optimal
size-connectance relationship for the assembly process.

Along the same line, one could look into the extension of our analysis beyond
the dynamics in Eq. (1), characterized by an interaction term depending
linearly on the $M$ matrix. This dynamics is largely employed in the
literature, but many authors also use its non-linear extension, commonly
termed Holling type II. The physical difference between the two frames lies in
the interaction (or handling) time: it is neglected in the linear case, but
not so in the other. Therefore, the magnitude of such a timescale introduces a
bound on the validity of the theory we presented. This is certainly available
when such a time can be approximated to zero. For finite non-zero values, one
may proceed perturbatively provided they are small. Still, it is possible to
show that as handling time tends the infinity, the resulting equation may be
cast in a linear form. Nevertheless, the effect of intermediate values of the
handling time represents an open question.

##  Appendix A Synthetic Network Generation

We start the analysis of the trade-offs between stability and feasibility for
a big ensemble of synthetic networks with the three patterns of interest, with
varying network size and connectance. The number of rows and columns nodes was
chosen within the range $[30,250]$, with $A=P$, and the connectance was mapped
onto the range $[0.04,0.4]$ such that the resulting networks were comparable
with empirical ecological networks.

We generated three different sets of synthetic networks, each one with purely
nested, modular, and in-block nested patterns, respectively. For this purpose,
we have employed a network generation model, introduced in Solé-Ribalta  _et
al._ (2018), to generate modular and in-block nested synthetic networks,
having a fixed block size ($N_{B_{A}}=N_{B_{P}}=15$ nodes) and an increasing
number of blocks, i.e., increasing network size. In the case of the nested
matrices, we employed the network generation model by keeping fixed the number
of blocks ($B=1$) and increasing its size $N_{B_{A}}=N_{B_{P}}$.

Besides the number of communities $B$, and the number of rows and columns
nodes, the model relies on three additional parameters: noise regarding the
existence of interactions outside species communities $\mu\in[0,1]$, noise
regarding interactions outside a perfectly nested structure $p\in[0,1]$ and a
shape parameter that controls the slimness of the nested structure
$\xi\in[1,\infty]$. These parameters were varied accordingly for each set of
matrices, e.g., for the generation of purely modular networks, we fixed $p=1$
and smoothly modified the rest of the parameters to obtain networks covering a
wide range of modularity values. The same scheme was applied for the
generation of purely nested, in-block nested networks, respectively.

Last but not least, given though $\xi$ affects the overall network connectance
(total number of existing species interactions), in some cases, we randomly
removed or added links to the networks until the desired connectance values
were met. The other parameters $p$ and $\mu$ did not alter the density of the
network. At the end of the process, we generated around $2\times 10^{5}$
networks for each one of the three patterns of interest.

##  Appendix B Nested, Modularity and In-Block nestedness measures

The present appendix is devoted to specify the measures we employ to quantify
the emergence of structural pattern in mutualistic networks, for instance to
build the simplex in Fig. 5.

### Nestedness

The notion of nestedness was coined over a decade ago in biology to
characterize the spatial distribution of biotas in isolated, fragmented
habitats Atmar and Patterson (1993). A perfect nested configuration (see Fig.
1, left)is observed when specialists (nodes with low connectivity) interact
with proper subsets of those species interacting with generalists (nodes with
high connectivity). To quantify this pattern we recall the global nestedness
fitness $\mathcal{N}$ introduced by Solé-Ribalta et al. in Solé-Ribalta  _et
al._ (2018), which is defined as

$\displaystyle\mathcal{N}=$
$\displaystyle\frac{2}{N_{A}+N_{P}}\sum^{N_{A}}_{i,j}\left[\frac{O_{i,j}-\langle
O_{i,j}\rangle}{k_{j}(N_{A}-1)}\Theta(k_{i}-k_{j})\right]$ $\displaystyle+$
$\displaystyle\frac{2}{N_{A}+N_{P}}\sum^{N_{P}}_{l,m}\left[\frac{O_{l,m}-\langle
O_{l,m}\rangle}{k_{m}(N_{P}-1)}\Theta(k_{l}-k_{m})\right],$ (11)

where $O_{i,j}$ (or $O_{l,m}$) measures the degree of links overlap between
rows node $i$ and $j$ (or coloumns $l$ and $m$) node pairs; the quantity
$k_{i}$ corresponds to the degree of the $i$ node; $\Theta(\cdot)$ is a
Heaviside step function that guarantees that we only compute the overlap
between pair of nodes when $k_{i}\geq k_{j}$. Finally, $\langle
O_{i,j}\rangle$ represents the expected number of links between row nodes $i$
and $j$ in the null model, and is equal to $\langle
O_{i,j}\rangle=\frac{k_{i}k_{j}}{N_{A}}$. This measure is in the tradition of
other overlap measures, such as NODF Almeida-Neto  _et al._ (2008); Ulrich
_et al._ (2009).

### Modularity

A modular networks configuration (Fig. 1, center) implies the existence of
well-connected subgroups, which may be identified given the right heuristics
to do so. Here, we search a (sub)optimal modular partition of the nodes
through a community detection analysis Barber (2007); Fortunato (2010). For
this goal, we apply the extremal optimisation algorithm Duch and Arenas (2005)
(along with a Kernighan-Lin refinement procedure Kernighan and Lin (1970)) to
maximise Barber’s modularity $\mathcal{Q}$ Barber (2007),

$\mathcal{Q}=\frac{1}{L}\sum^{N_{A}}_{i=1}\sum^{N_{A}+N_{P}}_{j=N_{A}+1}\left(\tilde{a}_{ij}-\tilde{p}_{ij}\right)\delta(\alpha_{i},\alpha_{j}),$
(12)

where $L$ is the number of interactions (links) in the network,
$\tilde{a}_{ij}$ denotes the existence of a link between nodes $i$ and $j$,
$\tilde{p}_{ij}=k_{i}k_{j}/L$ is the probability that a link exists by chance,
and $\delta(\alpha_{i},\alpha_{j})$ is the Kronecker delta function, which
takes the value $1$ if nodes $i$ and $j$ are in the same community, and $0$
otherwise.

### In-Block Nestedness

In-block nestedness configuration refers to a mesoscale hybrid pattern, in
which the network presents a modular structure but the interactions within
each module are nested (Fig. 1, right). In order to compute the amount of in-
block nestedness characterizing a given network we recall the objective
function introduced in Solé-Ribalta  _et al._ (2018), which employ the same
optimization algorithms used to maximize modularity. The in-block nestedness
objective function may be written as

$\displaystyle\mathcal{I}=$
$\displaystyle\frac{2}{N_{A}+N_{P}}\sum^{N_{A}}_{i,j}\left[\frac{O_{i,j}-\langle
O_{i,j}\rangle}{k_{j}(C_{i}-1)}\Theta(k_{i}-k_{j})\delta(\alpha_{i},\alpha_{j})\right]$
$\displaystyle+$
$\displaystyle\frac{2}{N_{A}+N_{P}}\sum^{N_{P}}_{l,m}\left[\frac{O_{l,m}-\langle
O_{l,m}\rangle}{k_{m}(C_{l}-1)}\Theta(k_{l}-k_{m})\delta(\alpha_{l},\alpha_{m})\right],$
(13)

where $C_{i}$ represents the size of the $i$ block. Note that, by definition,
$\mathcal{I}$ reduces to $\mathcal{N}$ when the number of blocks is $1$. This
explains why the right half of the ternary plot (Fig. 5) is necessarily empty:
$\mathcal{I}\geq\mathcal{N}$ implies that $f_{\mathcal{I}}\geq
f_{\mathcal{N}}$. On the other hand, an in-block nested structure exhibits
necessarily some level of modularity, but not the other way around. This
explains why the lower-left area of the simplex is empty as well (see Palazzi
_et al._ (2019a) for details).

## References

  * Atmar and Patterson (1993) W. Atmar and B. D. Patterson, Oecologia 96, 373 (1993).
  * May (1972) R. May, Nature 238, 413 (1972).
  * Bascompte  _et al._ (2003) J. Bascompte, P. Jordano, C. J. Melián, and J. M. Olesen, Proceedings of the National Academy of Sciences 100, 9383 (2003).
  * Mariani  _et al._ (2019) M. S. Mariani, Z.-M. Ren, J. Bascompte, and C. J. Tessone, Physics Reports 813, 1 (2019), nestedness in complex networks: Observation, emergence, and implications.
  * Olesen  _et al._ (2007) J. M. Olesen, J. Bascompte, Y. L. Dupont, and P. Jordano, Proceedings of the National Academy of Sciences 104, 19891 (2007).
  * Newman (2006) M. E. Newman, Proceedings of the national academy of sciences 103, 8577 (2006).
  * Lewinsohn  _et al._ (2006) T. M. Lewinsohn, P. Inácio Prado, P. Jordano, J. Bascompte, and J. M. Olesen, Oikos 113, 174 (2006).
  * Okuyama and Holland (2008) T. Okuyama and J. N. Holland, Ecology Letters 11, 208 (2008).
  * Thébault and Fontaine (2010) E. Thébault and C. Fontaine, Science 329, 853 (2010).
  * Allesina and Tang (2012) S. Allesina and S. Tang, Nature 483, 205 (2012).
  * Stouffer and Bascompte (2011) D. B. Stouffer and J. Bascompte, Proceedings of the National Academy of Sciences 108, 3648 (2011).
  * Staniczenko  _et al._ (2013) P. P. Staniczenko, J. C. Kopp, and S. Allesina, Nature communications 4, 1391 (2013).
  * Allesina  _et al._ (2015) S. Allesina, J. Grilli, G. Barabás, S. Tang, J. Aljadeff, and A. Maritan, Nature Communications 6, 7842 (2015).
  * Grilli  _et al._ (2016) J. Grilli, T. Rogers, and S. Allesina, Nature Communications 7, 12031 (2016).
  * Bastolla  _et al._ (2009) U. Bastolla, M. A. Fortuna, A. Pascual-García, A. Ferrera, B. Luque, and J. Bascompte, Nature 458, 1018 (2009).
  * Saavedra  _et al._ (2015) S. Saavedra, R. Rohr, J. Olesen, and J. Bascompte, Ecology and evolution 6, 1007 (2015).
  * Suweis  _et al._ (2013) S. Suweis, F. Simini, J. R. Banavar, and A. Maritan, Nature 500, 449 (2013).
  * Tilman and Downing (1994) D. Tilman and J. Downing, Nature 367, 363 (1994), copyright: Copyright 2018 Elsevier B.V., All rights reserved.
  * McCann (2000) K. S. McCann, Nature 405, 228 (2000).
  * Palazzi  _et al._ (2019a) M. Palazzi, J. Borge-Holthoeffer, C. Tessone, and A. Solé, J.R. Soc. Interface 16, 20190553 (2019a).
  * Gracia-Lázaro  _et al._ (2018) C. Gracia-Lázaro, L. Hernández, J. Borge-Holthoefer, and Y. Moreno, Scientific reports 8, 1 (2018).
  * Stella  _et al._ (2018) M. Stella, S. Selakovic, A. Antonioni, and C. S. Andreazzi, Elife 7, e32814 (2018).
  * Pilosof  _et al._ (2017) S. Pilosof, M. A. Porter, M. Pascual, and S. Kéfi, Nature Ecology & Evolution 1, 1 (2017).
  * Kéfi  _et al._ (2016) S. Kéfi, V. Miele, E. A. Wieters, S. A. Navarrete, and E. L. Berlow, PLoS biology 14, e1002527 (2016).
  * García-Callejas  _et al._ (2018) D. García-Callejas, R. Molowny-Horas, and M. B. Araújo, Oikos 127, 5 (2018).
  * Hutchinson  _et al._ (2019) M. C. Hutchinson, B. Bramon Mora, S. Pilosof, A. K. Barner, S. Kéfi, E. Thébault, P. Jordano, and D. B. Stouffer, Functional Ecology 33, 206 (2019).
  * Fortuna  _et al._ (2010) M. A. Fortuna, D. B. Stouffer, J. M. Olesen, P. Jordano, D. Mouillot, B. R. Krasnov, R. Poulin, and J. Bascompte, Journal of Animal Ecology 79, 811 (2010).
  * Flores  _et al._ (2011) C. O. Flores, J. R. Meyer, S. Valverde, L. Farr, and J. S. Weitz, Proceedings of the National Academy of Sciences 108, E288 (2011).
  * Flores  _et al._ (2013) C. O. Flores, S. Valverde, and J. S. Weitz, The ISME Journal 7, 520 (2013).
  * Beckett and Williams (2013) S. J. Beckett and H. T. Williams, Interface Focus 3, 20130033 (2013).
  * Solé-Ribalta  _et al._ (2018) A. Solé-Ribalta, C. J. Tessone, M. S. Mariani, and J. Borge-Holthoefer, Physical Review E 96, 062302 (2018).
  * Palazzi  _et al._ (2019b) M. J. Palazzi, J. Cabot, J. L. Cánovas Izquierdo, A. Solé-Ribalta, and J. Borge-Holthoefer, Scientific Reports 9, 2045 (2019b).
  * Cai  _et al._ (2020) W. Cai, J. Snyder, A. Hastings, and R. M. D’Souza, Nature Communications 11 (2020), https://doi.org/10.1038/s41467-020-19154-5.
  * Pinheiro  _et al._ (2019) R. B. Pinheiro, G. M. Felix, C. F. Dormann, and M. A. Mello, Ecology 100, e02796 (2019).
  * Gibbs  _et al._ (2018) T. Gibbs, J. Grilli, T. Rogers, and S. Allesina, Phys. Rev. E 98, 022410 (2018).
  * Svirezhev and Logofet (1983) I. Svirezhev and D. O. Logofet,  _Stability of biological communities_ (Mir Publisher, 1983).
  * Grilli  _et al._ (2017) J. Grilli, M. Adorisio, S. Suweis, G. Barabás, J. R. Banavar, S. Allesina, and A. Maritan, Nature Communications 8, 14389 (2017).
  * Ribando (2006) J. M. Ribando, Discrete and Computational Geometry 36, 479 (2006).
  * Genz and Bretz (2009) A. Genz and F. Bretz,  _Computation of Multivariate Normal and t Probabilities_ (Springer-Verlag Berlin Heidelberg, 2009).
  * (40) “Web of life: ecological networks database,” http://www.web-of-life.es/.
  * Almeida-Neto  _et al._ (2008) M. Almeida-Neto, P. Guimarães, P. R. Guimarães, R. D. Loyola, and W. Ulrich, Oikos 117, 1227 (2008).
  * Ulrich  _et al._ (2009) W. Ulrich, M. Almeida-Neto, and N. J. Gotelli, Oikos 118, 3 (2009).
  * Barber (2007) M. J. Barber, Physical Review E 76, 066102 (2007).
  * Fortunato (2010) S. Fortunato, Physics Reports 486, 75 (2010).
  * Duch and Arenas (2005) J. Duch and A. Arenas, Physical review E 72, 027104 (2005).
  * Kernighan and Lin (1970) B. W. Kernighan and S. Lin, The Bell system technical journal 49, 291 (1970).