# Roles of Dynamic State Estimation in Power System Modeling, Monitoring and
Operation

IEEE Task Force on Power System Dynamic State and Parameter Estimation
Junbo Zhao (TF Chair), , Marcos Netto, , Zhenyu Huang, , Samson Shenglong Yu,
, Antonio Gómez-Expósito, , Shaobu Wang, , Innocent Kamwa, , Shahrokh
Akhlaghi, , Lamine Mili, , Vladimir Terzija, , A. P. Sakis Meliopoulos, ,
Bikash Pal, , Abhinav Kumar Singh, , Ali Abur, , Tianshu Bi, , Alireza Rouhani

###### Abstract

Power system dynamic state estimation (DSE) remains an active research area.
This is driven by the absence of accurate models, the increasing availability
of fast-sampled, time-synchronized measurements, and the advances in the
capability, scalability, and affordability of computing and communications.
This paper discusses the advantages of DSE as compared to static state
estimation, and the implementation differences between the two, including the
measurement configuration, modeling framework and support software features.
The important roles of DSE are discussed from modeling, monitoring and
operation aspects for today’s synchronous machine dominated systems and the
future power electronics-interfaced generation systems. Several examples are
presented to demonstrate the benefits of DSE on enhancing the operational
robustness and resilience of 21st century power system through time critical
applications. Future research directions are identified and discussed, paving
the way for developing the next generation of energy management systems.

###### Index Terms:

Dynamic state estimation, Kalman filtering, low inertia, monitoring, parameter
estimation, power system stability, synchronous machines, synchrophasor
measurements, converter interfaced generation, static state estimation.

##  I Introduction

Dynamic state estimation (DSE) [1] is going to be very useful for time
critical monitoring, control, and protection of future electric power grids.
This is largely due to the changes in generation mixes and load compositions,
particularly the increasing penetration of intermittent, stochastic and power
electronics-interfaced non-synchronous renewable generation and distributed
energy resources (DERs) [2]. In this context, where the system operating point
changes more often and more rapidly, tracking the system dynamic state
variables is of critical importance. Recall that dynamic state variables are
the ones associated with the time derivatives in the set of differential-
algebraic equations describing power system dynamics [3]. The dynamic state
variables of the synchronous generators correspond to their electromagnetic
and electromechanical processes, as well as their controllers [4]. As for the
non-synchronous generations, the dynamic state variables are associated with
the primary source of energy, e.g. solar photovoltaic arrays, batteries and
wind turbines, and their controllers. Note that the phenomena driving the
primary source of energy in this case are diverse. Examples include frequency
at the point of common coupling between the power converter and the electric
grid, the electric current in the converter, and the pitch angle of wind
turbines. Furthermore, there exists an important challenge in the development
and maintenance of accurate models for power electronics-interfaced devices.
DSE can be developed to validate the models and to estimate unknown or
incorrect parameters.

Despite being initially mentioned in the 1970s [5], it was only in recent
years that the power system community has picked up the momentum in DSE
research. Part of the reason was the lack of appropriate metering
infrastructure, like phasor measurement units (PMUs) and merging unit (MUs)
that are being widely deployed to capture the appropriate dynamics in power
systems [6, 7]. Feasibility studies using PMU measurements for DSE are
reported in [8], and subsequently, various Kalman filtering techniques, such
as Extended Kalman Filter (EKF), Unscented Kalman Filter (UKF) [9], Ensemble
Kalman Filter (EnKF), Particle Filter (PF) and their variants [10] have been
applied to DSE. Data-driven DSE [11], observability analysis to guide
measurement selection [12] and the enhancement of robustness against bad data
and parameter errors are also developed in [13, 14, 15]. Readers may refer to
[4] for a comprehensive summary of DSE algorithms.

Several online tasks can benefit from DSE but have not adequately deliberated
the need for it. An example is the dynamic security assessment (DSA). Today’s
DSA tool assumes a steady-state initial condition [16], which would yield
inaccurate results in the presence of more fluctuations from DERs and loads.
While DSE could directly estimate the dynamic states and therefore provide
more accurate initialization conditions [17]. Furthermore, for synchronous
generators, numerous developed control schemes rely on power system
stabilizers (PSS) to damp out oscillations. These control schemes use the
rotor’s speed as input. The direct use of measurements obtained from a meter
installed on the shaft of the machine is not reliable; instead, the current
practice is to use the compensated frequency [18], which is calculated using
the voltage and current measured at the generators’ terminal. However, the
calculation of the compensated frequency is significantly affected under
transient conditions [19], which leads to sub-optimal control response of the
PSS. If DSE is deployed for frequency estimation, this issue can be
effectively resolved. The rate of change of frequency (ROCoF) estimation can
also be formulated in the state- space model for DSE [20]. Note that, in these
examples, traditional static state estimation (SSE) could not provide reliable
system dynamic states and quantities. Therefore, it is of critical importance
to understand the role of DSE in power systems from various aspects. While SSE
has been in energy management systems (EMS) for decades and its roles and
implementation requirements are well established, the same does not apply for
DSE. There is a large gap in understanding different roles of SSE and DSE and
how DSE can be implemented in practice.

This paper summarizes the joint efforts of the Task Force in clarifying the
roles of DSE from the power system modeling, monitoring, and operation
perspectives. Comparisons between SSE and DSE in practical implementations are
thoroughly discussed, including measurement, model requirements, software
support and potential applications based on the estimation results.
Illustrative examples and summary frameworks have been presented to appreciate
the roles of DSE. The future applications are discussed to shed light on the
transition from today’s EMS to its next generation.

The remainder of this paper is organized as follows. Section II compares DSE
and SSE from the practical implementation requirement and application aspects.
Section III summarizes the role of DSE in different applications. Conclusions
and directions for future research are given in Section IV.

TABLE I: Comparison Between SSE and DSE. | SSE | DSE
---|---|---
Measurements | From the SCADA | From PMU/DFR
|  (every $\sim$2-10 sec.) |  (every 1/30$\sim$1/240 sec.)
Observability |  Binary111Topological observability analysis provides a binary answer as the first byproduct, but it gives additional valuable information, such as which buses (islands) are observable or which pseudo-measurements should minimally be added to restore full observability. Numerical observability analysis, being based on the factorization of the Jacobian or Gain matrices, provides a “spectrum” or range of observability answers, depending on the condition number of the matrix being factorized. A network can be topologically observable but not algebraically (numerically), owing to abnormal network parameters or measurement weights. In some cases, a network can be numerically observable for a certain combination of weights, but not for others. This is not exactly “binary”. | Time-varying
| (observable or not) | (Strong/weak/not observable)
Update speed | 1 snapshot | 1 prediction + 1 filtering
|  (every $\sim$2-10 sec.) | step of the Kalman filter
|  |  (every $<$ 1/30 sec.)
Models | Algebraic power | Differential-algebraic
| flow equations | equations
Framework | Mostly centralized | Centralized &
| or distributed | distributed/decentralized
Outputs | Algebraic variables | Dynamic variables
| (voltage magnitudes | (machine/ dynamic load/
| and angles) | DERs dynamic variables)
Applications | Monitoring and control | Monitoring (operator in
| (operator in the loop) | the loop), control and
|  | adaptive protection
|  | (operator out of the loop)

##  II Comparative Overview of SSE and DSE

SSE has become a widely used tool in today’s EMS, while DSE is a new tool for
the industry and system operators. It is essential to clarify their
implementation and functionality differences and, at the same time, enable a
clear path transition from SSE-based EMS to the future DSE-based EMS with
power electronics-dominated power systems.

###  II-A Implementation Differences

SSE and DSE have different requirements with respect to measurements, models,
observability, execution rate, outputs and applications. These differences are
summarized in Table I. SSE mostly relies on SCADA measurements that are
updated every 2-5s and some PMU measurements to gain more redundancy. But
since SCADA measurements are not synchronized while those of PMUs are,
effectively integrating those two sources of data needs to be taken care of.
By contrast, for DSE, fast and synchronized measurements with a reporting rate
of 30 to 240 samples per second, are used and they might come from PMUs and
digital fault recorders (DFRs). Furthermore, there is a significant difference
in the observability theory for SSE and DSE [12, 21]. For SSE, the
topological- or numerical-based observability analysis typically determines
whether the system is observable or not. If the system is unobservable,
observable islands may be determined [22]. Furthermore, in the presence of
ampere measurements, where multiple solutions may exist depending on the
nature of the remaining measurements and the loading point, the observability
analysis can also inform whether there are multiple solutions or not to the
SSE problem. A network can be uniquely observable for heavily loaded cases,
but “unobservable” (i.e., undefined Jacobian) when the ampere measurement
value is nearly zero (unloaded line), which leads to 0/0 indetermination.
Thus, the answer depends on the line loading but is typically binary. By
contrast, in DSE, one may refer to strongly or weakly or not observable
systems. One way to quantify this is to compute the smallest singular value of
the observability matrix from the Lie-derivatives. Higher (lower) values of
the smallest singular value of the observability matrix indicates stronger
(weaker) observability for a given measurement set [12]. Furthermore, due to
the nonlinear and time-dependent nature of the problem, observability results
are time-varying. It is worth noting that under certain conditions, the system
might be unobservable, but still detectable [23] for DSE. Detectability is a
slightly weaker notion than observability. A system is detectable if all the
unobservable states are stable. Also, each prediction-correction step of the
DSE must be numerically solved faster than the PMUs/DFRs/MUs scan rate,
thereby posing some challenges on the computational power. To address this
issue, decentralized/distributed DSE or parallel computing technique for
centralized DSE is usually suggested [4]. The outcome of SSE and DSE are also
different. SSE provides estimates of the bus voltage magnitudes and phase
angles, which are the algebraic variables; on the other hand, DSE provides
estimates of the dynamic state variables, such as those associated with
generators/dynamic loads/DERs. There is also joint DSE that estimates dynamic
and algebraic variables simultaneously [4]. For some PMU observable networks,
linear state estimation (LSE) that keeps up with the PMU refreshing rate is
developed [24]. However, it does not track the actual system dynamics.

When implementing SSE and DSE, the models used are different. In SSE, the
generators and loads are simply modeled by power injections and, hence, the
system model used for SSE is represented using algebraic equations. While for
DSE, the generators/dynamic loads/DERs, etc., and their controllers are
represented by a set of DAEs. Note that it is not required to have PMUs
installed at each generator terminal; if the generator terminal is observable
via a local LSE, the DSE can be implemented.

Figure 1: Functional schematic comparisons between SSE and DSE.

###  II-B Functionality Differences

From the functionality perspective, both SSE and DSE have two essential roles,
namely improving the targeted system visibility and validating the used data
and models, which are summarized in Fig. 1. Since SSE is the snapshot-based
static view of the system algebraic variables, its outputs allow the
calculations of power injections and power flows across the systems. This
provides the database for many applications in today’s EMS, including system
operations, such as power flow, optimal power flow, load forecasting, etc., as
well as security assessment, i.e., static contingency and static voltage
stability analysis. By contrast, DSE relies on fast sampled and time
synchronized measurements to track system dynamic changes, thus providing
dynamic visibility, such as dynamic states, anomaly detection, unknown inputs
estimation, etc. [25, 26, 27]. This dynamic information further allows other
related applications, i.e., DSA [17], rotor angle stability assessment [27],
adaptive power system protection [28] and so on. In terms of data and model
validation for SSE, it relies on the measurement residual-based statistical
test or robust estimation criteria for bad data processing and measurement
model parameter error identifications. The bad data are caused by gross
measurement errors or data transfer errors while model parameter errors can be
due to transmission line parameters and topology inaccuracies. Since DSE
relies on a set of DAEs, besides those parameter errors previously mentioned,
it also encounters in equations that govern the system dynamics, e.g.,
generator models and their related controllers. These errors are classified
into two categories: innovation outliers, which are caused by parameter and/or
input errors in the set of DAEs; and structural outliers, which are caused by
a model structure that does not reflect well the system dynamics [29]. Note
that the DSE outputs can be used to calculate the system algebraic variables
as well. Therefore, the SSE functions can benefit from DSE. And vice versa,
the information provided by the SSE can be used to warm start the DSE process
before an event occurs.

_Remark:_ with the increasing deployment of PMUs, many high voltage
transmission systems are already observable with only PMU measurements, such
as that of the Dominion Energy, the 765/345/230 kV power grid in New York, and
the 345 kV power grid in New England. As a result, the legacy nonlinear SSE
algorithms can now be reformulated as the LSE [24]. As long as the
computational power is sufficient, the LSE can be updated with the scan rate
of PMUs but its results are also restricted to algebraic variables-based
applications. Compared with the SCADA-based SSE, an LSE allows for faster
real-time contingency analysis and voltage stability assessment, besides area
angle limit monitoring. By comparing LSE with DSE, LSE does not touch the
dynamic equations and can not provide the real-time picture of system dynamic
states and the related controllers. For example, the rotor speed information
is not available from LSE and thus many applications based on them cannot be
done. The importance of dynamic states for power system monitoring, visibility
and operation will be highlighted in Section III. In summary, the fundamental
differences in terms of potential applications are that LSE only deals with
algebraic variables while DSE provides both dynamic and algebraic variables of
the system.

###  II-C Practical Implementation of DSE

There are several vendors who offer commercial software for SSE implementation
while there is no commercial software for DSE. However, we can leverage the
capabilities of existing commercial tools and enable the DSE implementation.
The practical implementation of DSE is shown in Fig. 2, where the system is
modeled by the following state-space model after time discretization of DAEs,
and considering the available measurements:

$\displaystyle\bm{x}_{k}$
$\displaystyle=\bm{f}(\bm{x}_{k-1},\bm{y}_{k-1},\bm{u}_{k},\bm{p}_{k})+\bm{w}_{k},\;\mathbb{E}\left[\bm{w}_{k}\bm{w}_{k}^{\top}\right]=\bm{Q}_{k}$
(1) $\displaystyle\bm{z}_{k}$
$\displaystyle=\bm{h}(\bm{x}_{k},\bm{y}_{k},\bm{u}_{k},\bm{p}_{k})+\bm{v}_{k},\;\mathbb{E}\left[\bm{v}_{k}\bm{v}_{k}^{\top}\right]=\bm{R}_{k}$
(2)

where $\bm{x}_{k}$ and $\bm{y}_{k}$ represent system dynamic and algebraic
state vectors, respectively; $\bm{z}_{k}$ is the measurement vector from PMUs
or DFR or MU; $\bm{u}_{k}$ is the system input vector that drives the state
transition; $\bm{p}_{k}$ denotes the system parameters; $\bm{f}$ and $\bm{h}$
are vector-valued nonlinear functions; $\bm{w}_{k}$ and $\bm{v}_{k}$ are the
system process and measurement error vector, respectively, with covariance
matrices $\bm{Q}_{k}$ and $\bm{R}_{k}$.

Figure 2: Implementation of DSE using existing commercial tools.

The implementation of DSE has two main steps, namely the state prediction via
(1) to obtain predicted state at time $k$ and the state correction that
integrates the state prediction and measurements in (2) for filtering. It is
interesting to note that the widely used commercial software, such as PSSE/
PSLF/ PowerWorld/ DSATools have the so-called event playback functions
available. This play back function allows us to implement equations (1)-(2)
automatically. The remaining work is to rely on Kalman filter framework and
integrate the state prediction and the measurements together for the
estimation of dynamic states. To fully leverage this playback function, it is
suggested to use nonlinear Kalman filter techniques that do not require
linearization, such as UKF, EnKF, PF and their variants. Another useful
implementation of DSE is to take unknown dynamic states or control inputs as
additional state variables for joint estimation. This is because there may be
lack of accurate controller information or the control models may have poor
quality. Therefore, it can be concluded that only the filtering step needs to
be coded along with the playback function for practical DSE implementation. It
should be noted that we do not need to wait for an event to trigger the DSE.
As long as the measured vector $\bm{u}_{k}$, typically the voltage magnitudes
and angles/frequency, is fed into the playback function, the DSE can be
implemented to continuously monitor dynamic states. Note that the
observability and required measurements should be provided when implementing
DSE [12].

##  III Roles of DSE in Power Systems Modeling, Monitoring and Operation

Three major application areas that can potentially benefit from DSE have been
identified and elaborated in this paper, namely: a) modeling, b) monitoring
via the enhanced dynamic visibility, and c) operation. Each of them is
discussed below.

###  III-A Modeling

Following the 2003 Northeastern U.S. and Eastern Canada blackout, the North
America Electric Reliability Corporation (NERC) has been developing standards
for periodical model validation [30]. In the existing commercial software, the
“event playback” function is leveraged to validate dynamic models using PMU
data. The key idea is to take the measured generator terminal voltage
magnitude/phase angle or frequency as model inputs to obtain model outputs of
real power $P$ and reactive power $Q$. These responses are compared with the
measured $P$ and $Q$ to validate the model adequacy [31, 32, 33]. When the
model of the subsystem is accurate, the simulation responses should match
statistically well the actual responses. By contrast, when significant
mismatches are observed, the model is considered inadequate and parameter
calibration is needed. This process can be found in the upper part of Fig. 3,
where $\bm{e}_{k}$ is the parameter error with covariance matrix $\bm{W}_{k}$
and $\bm{g}$ is the parameter regression function. The difference between
Figs. 2 and 3 is that additional parameters are augmented with the states for
joint estimation. Note that the criterion for model validation is a byproduct
of DSE, namely the innovation vector-based statistical test. Furthermore, the
targeted dynamic components can be synchronous machines, wind farms, power
electronics-interfaced DERs [34] that can be described by DAEs.

Figure 3: DSE-enabled model validation and calibration.

To calibrate the system models, both the system model structure and parameters
need to be carefully investigated. Traditionally, they can be found from
manufacturer databases. However, parameters might drift during operations due
to a variety of factors, e.g., environmental changes, the aging process, and
coupling effects. Taking generators offline for testing and parameter
calibration is costly and affects the system reliability. Therefore, online
calibration has been an important application of DSE. Since not all parameters
are identified using the present measurement sets, it is, therefore, critical
to select appropriate candidate parameters before starting a parameter
calibration procedure. The calibration typically includes four main steps [31,
32], namely 1) initial checks, 2) sensitivity analysis, 3) parameter
estimation and 4) parameter validation using other disturbances. The first
step is to eliminate obvious errors in parameters that are not realistic, such
as fractional values for the model flags that should be integers, parameter
values outside the normal ranges, swapped values for limiters, and incorrect
statuses of controllers. The sensitivity analysis provides a guideline in
determining parameter sensitivity for different disturbances and identifying
the candidate parameters for calibration. The selected parameter vector
$\bm{p}$ is augmented with the original state vector for joint estimation as
shown in the lower part of Fig. 3. The joint DSE methods are based on, e.g.,
EKF, UKF, constrained UKF, and EnKF [35, 36, 37, 38]. After calibration, it is
necessary to perform model validation using different events so that the
identified parameters are not local optima. If the model deficiency for other
events is detected, the calibration process should be repeated until there are
no inconsistent responses between models and measurements. Besides the model
validation and calibration, the power system coherency identification would
also benefit from the estimated dynamic states, especially the rotor angles
and speeds [39, 40]. The identified coherency allows the development of model
reduction for stability assessment and islanding control. Another potential
application of DSE is the dynamic reduction of large power systems [41], where
the equivalent model parameters can be estimated together with the machine
dynamic states. The dynamic load parameter identification and the unknown
machine parameters, such as inertia, sub-transient reactance, governor dead-
band, etc., can also be cast into the state-space formula. Then, various DSE
algorithms can be utilized for unknown parameter estimation [42, 43, 44].

###  III-B Monitoring

One of the key DSE applications is to provide enhanced dynamic visibility for
power system monitoring. Part of the dynamic visibility of DSE has been shown
in Fig. 1 but a more detailed summary includes seven key functions: i) dynamic
state trajectory tracking, ii) oscillation monitoring, iii) bus frequency,
ROCOF and center of inertia (COI) frequency estimation, iv) data quality
detection and correction, including cyber attacks, v) unknown control inputs
identification, vi) anomaly detection and vii) other applications that involve
with larger time constants than the electromechanical dynamics. The dynamic
state trajectory tracking is a natural product of DSE as it provides the time
series of the generator and its controller states in the presence of system
disturbance. Among these, the generator rotor speed and angle have been widely
used in power system, such as oscillation detection [45]. For example, in
damping-torque-based forced oscillation source location, the rotor speed, and
angle are required [46]. The estimated machine rotor speeds also allow us to
resort to the frequency divider for bus frequency estimation [47]. This
interesting result reveals that if the system machine rotor speeds are
available, all bus frequencies can be estimated. Since the number of machines
is much smaller than the number of buses, it significantly reduces the
requirement of PMU installations for bus frequency monitoring. The real-time
COI frequency plays an important role in power system stability analysis and
control. Via the DSE outputs, we can obtain the COI frequency at the same
refreshing rate of the PMU measurements. This provides the online reference
frequency for control [48]. Note that a widely used industry practice is to
leverage generator terminal voltage and current measurements to calculate an
approximated internal machine rotor speed. However, its accuracy is
questionable as can be seen in the following example: a three-phase short-
circuit occurs at Bus 16 of the IEEE 39-bus system on $t$=0.5s and the fault
is cleared after one cycle by opening the transmission line connecting buses
16 and 17. The machines are modeled using the subtransient model with a DC1A
exciter, PSS1A stabilizer, and IEEE-G3 governor. Fig. 4 shows the estimated
generator 5 internal rotor speed by DSE, numerical derivative with low pass
filter and numerical derivative with washout filter [49] of the derived
generator angle from the terminal voltage and current measurements. As
expected, the numerical derivative of the approximated rotor angle suffers
from numerical errors at the moment the event happens; numerical derivative
with washout filter is able to address that but at the cost of modifying the
time-domain response. By contrast, the DSE can accurately track the rotor
speed. On the other hand, ROCoF has been utilized as an important parameter
for system protection, especially for the inverter-based generations. By
treating it as a dynamic state instead of an algebraic variable, an accurate
estimation of it can be achieved [20].

Figure 4: Comparisons of different methods for internal generator frequency
estimation.

The measurements are always subject to noise and even large errors caused by
communications, instrument transformer saturation, cyber-attacks, etc. It is
important that DSE naturally filters out the measurement noise and provides
more accurate data for other applications. Since the field measurements
typically follow non-Gaussian distributions [50], robust DSE developed based
on the robust statistics theory is needed [51]. While for randomly occurred
bad measurements, both normalized residual statistical tests and robust
detection are developed. However, it is shown that the threshold for the
traditional normalized residual statistical test is not analytical and depends
on different systems characteristics. Also, it does not provide reliable
outputs in the presence of unknown measurement noise statistics. This is not
the case for the robust estimation that automatically detects and suppresses
bad data. For the cyber-attack scenario, robust DSE with high breakdown points
or novel machine learning-aided DSE might offer a solution [52]. This is still
an open area that requires further investigation.

Another important visibility is on the excitation system that could
significantly affect the system stability. Although PMU measurements might be
used to calculate the excitation voltage in certain cases [53], this cannot be
generalized to the modern brushless excitation systems. Furthermore, under
stressed conditions, the excitation voltage is restricted by timer-based over-
excitation limiters, which may dramatically affect the system stability margin
[54]. Traditional DSA model typically does not capture if the state variables
corresponding to exciter output are saturated or not. When saturation occurs,
it is no longer a state variable and thus should be kept as an unknown input
estimation. To achieve this goal, the excitation voltage is taken as the
unknown variable for a two-stage estimation, where the estimation of input is
done first and then substituted into the original model for DSE. Both EKF [55,
56] and UKF [57] considering unknown excitation voltage are used for that and
the UKF has a better capability of handling the model nonlinearities and
saturation effect. These two approaches require the local generator frequency
to ensure the observability of unknown inputs. However, the local generator
frequency is typically not measured directly by the PMUs. To deal with that, a
correlation aided robust DSE for unknown input and state estimation is
proposed in [58]. Compared with the previous works, no generator frequency
measurement is required, and it has better robustness in dealing with bad
data. It is worth noting that the unknown (i.e., not measured) mechanical
torque of the turbine-governor system is also estimated using the approach in
[58].

Besides the unknown control inputs, there are also anomalies, such as
controller failures or malfunctions that can affect the system stability. For
example, if the over-excitation occurs, a voltage security issue may arise.
Furthermore, if the failure of the excitation system happens but without being
detected, the accuracy of the differential-algebraic equations will be
negatively impacted, leading to incorrect conclusions about system stability.
It is very important to detect these failures timely to avoid the risk of
exciter and voltage regulator damages. By relying on the estimated dynamic
states and checking the consistency between the control model outputs and
expected outputs, the over-excitation and abnormal mechanical power changes
can be detected [59]. The multiple models based DSE technique is also
developed to detect excitation failures in [25]. Those techniques may be
further extended to deal with the anomalies of other controllers, such as
governor, and power system stabilizer.

Besides monitoring electromechanical dynamics, there are other DSE-based
monitoring applications that use DAEs but with much higher time constants.
This is the case, for instance, of the increasingly important real-time
thermal rating of lines and cables, which is drastically changing the
customary way in which static security assessment is performed. Indeed, the
existing assets are reaching their conservative ampacity limits, typically
defined on a seasonal basis, without due consideration to meteorological
conditions (temperature, wind speed, etc.). The simultaneous, real-time
estimation of the external parameters arising in the thermal models of lines
and cables can help us improve the operation of underground cables [60, 61].

Figure 5: Influence of initial conditions on DSA results.

###  III-C Operation

SSE plays an important role in today’s power system operations, especially the
essential data base for contingency analysis, static voltage stability and
optimal power flow. Compared with SSE, DSE outputs would also offer new
benefits to today’s EMS. Two main categories are elaborated here, namely DSA
and dynamic stability assessment.

Model-based DSA relies on the initialization of dynamic models. The inaccurate
initialization of variables might lead to different conclusions about whether
the system is stable or not. In the legacy power systems, the operating point
evolves slowly in time and the SSE results can be used to calculate the
system’s initial conditions that are required by DSA tools [16]. In other
words, if the system is in equilibrium at 60Hz (or 50Hz), SSE is adequate for
initializing the differential equations, which is how DSA is done in today’s
EMS. When the system is not at equilibrium or not at 60Hz (or 50Hz), SSE can’t
be used to initialize the differential equations. Indeed, in future power
systems with high penetration of DERs, the variables will change more often
and more rapidly, rendering SSE too slow to provide DSA tools with an accurate
picture of the system states. But it is during this time period that DSA would
be even more important for power system operation. One example is the 2003
Northeast blackout, where the system frequency deviated from 60Hz for hours
and DSA would help to guide operator’s restorative actions. Model
initialization for DSA is therefore the first important role of DSE. This is
demonstrated through a simple numerical simulation carried out on the IEEE
39-bus system. The system is assumed to be operated under the scenario, where
the loads and DERs have large stochastic behaviors, yielding some oscillations
of the synchronous generators. The operator would like to assess if the system
is able to withstand a contingency. To this end, at $t$=1.3s, a three-phase
short-circuit is applied at Bus 28 and cleared after 30ms by opening the
transmission line connecting Buses 28 and 29. Two cases are considered: 1) the
model is initialized by using the estimated bus voltage magnitudes and angles
obtained from SSE at $t$=1.3s; 2) the model is initialized by using the
dynamic states obtained from DSE at $t$=1.3s. DSA is then performed in both
cases. The rotor angle of Generator 9 with respect to that of Generator 10 is
displayed in Fig. 5. It can be observed that the SSE-based DSA indicates that
the system loses stability while the actual system remains stable. By
contrast, the DSE-based DSA reflects true system behavior. The key insight is
that in the SSE-based initialization, the derivatives of dynamic state
variables are set to zeros based on the quasi-steady-state assumption, but the
actual system already deviates from a steady state due to the variations of
relative large loads and DERs. As a result, the calculated states using SSE-
based DSA are not accurate. This is not the case for DSE-based DSA as their
results contain the non-zero derivatives and can be directly used for non
quasi-steady-state initialization.

Figure 6: DSE for unstable system with high non-linearity.

Dynamic stability assessment typically relies on a database that is built from
offline simulations, wherein the evolution of the system dynamic states
following different events is recorded. Then, an online classifier algorithm,
e.g. decision tree or random forest, continuously compares real-time
measurements with the database. If an evolving instability is detected,
remedial control actions are initiated. This control philosophy is deployed in
many power systems around the world. Due to its dependence on offline
simulations, it is challenging to consider all possible operating scenarios
during offline simulations [62]. This is particularly important with the
increasing deployment of power electronics-interfaced DERs. DSE offers an
alternative to these databases as it can provide a picture of the system
dynamics in almost real-time, and can be used to compute transient stability
indices, e.g., rotor angle stability [63]. This alternative control philosophy
based on real-time information is yet in its infancy but seems to be
promising. Rotor angle stability can also be assessed via the sign of the
system’s maximal Lyapunov exponent [27], where the generator rotor angles
provided by DSE would be of importance to predict system stability. To further
demonstrate the capability of DSE in dealing with system instability, the same
three-phase fault for the results in Fig. 4 is used but with 30 cycles of
clearing time. This leads to system instability as shown in Fig. 6. The robust
UKF-based DSE [14] is utilized to estimate and track generator rotor angle and
speed. We observe that even on the verge of losing stability, where the system
response is highly nonlinear, DSE continues to track the evolution of the
system states that can be used for emergency/ in extremis control schemes. On
the other hand, the status and/or states of dynamic devices are playing a key
role in dynamic voltage stability, such as overexcitation limiters, and
dynamic loads. DSE can be appropriately formulated to monitor them and provide
insightful information on voltage stability conditions of the system. For
future power electronics-dominated systems, the frequency issue is of critical
concern due to the reduced system inertia. DSE can provide online frequency
information to enable more effective system frequency stability assessment.
However, there are very limited researches along these directions and more
investigations are needed.

###  III-D DSE for Power Electronics-Interfaced Renewable Generation
Visibility

Renewable energy sources (RESs) are typically integrated with the grid through
power electronics converters. To achieve reliable and cost-effective RESs
integration, dynamic visibility is needed. Compared with the conventional
synchronous generators, RESs exhibit a higher level of variability and
uncertainty. Furthermore, the power converters and the control loops of RESs
are very different from the synchronous machines. The power electronics have
much smaller time constants than those of the synchronous generators, yielding
much faster dynamics. As a result, the measurement system used for the
monitoring of RESs should have a faster sampling speed and broader bandwidths
for data transmission. It should be noted that besides the traditional voltage
stability, rotor angle stability and frequency stability, there is also
converter-induced stability [64]. For monitoring and control such fast
dynamics of the converters, DSE would be a good choice. For example, the UKF-
based DSE is developed to extract the fundamental components of the point of
common coupling voltage and load current for grid synchronization of PV
systems in [65]. The doubly fed induction generator (DFIG) wind generator’s
model can be put into the state-space form with their dynamic states estimated
by DSE to achieve visibility. The DSE can be implemented via UKF, particle
filter, and unscented particle filter [66, 67, 68]. Since there are errors
with the wind speed measurements, DSE for DFIG with unknown wind speeds is
addressed in [69]. The joint estimation of dynamic states and parameters of
permanent-magnet synchronous motor-based wind generator is investigated in
[70]. It should be noted that the relation of estimated states via DSE with
system dynamics and stability phenomena is still the subject of ongoing
research. Except for the visibility of dynamic states, it is also important to
detect any anomalies inside the RESs, such as control failures, erroneous
tripping actions, etc. This would provide timely information for operators to
take proactive controls for maintaining system stability. For example, due to
the erroneous frequency and ROCOF measurements provided by the phase-locked-
loop, several large solar farms have been incorrectly tripped in the 2016
California blue cut fire event [71]. The DSE would allow us to provide more
accurate frequency and ROCOF estimates and visibility of abnormal behaviors by
checking the model and measurement consistency. This is an open problem that
needs more investigations.

##  IV Conclusions and Future Work

This paper provides a comprehensive summary of the roles of DSE for power
system modeling, monitoring and operation. The relationships between SSE and
DSE have been identified from the implementation requirements and related EMS
functionalities. Several representative examples have been presented to
highlight the critical importance of DSE for future power systems, especially
for the development of the next-generation EMS.

Future research on DSE can be generally categorized into four key aspects:
data infrastructure for DSE, DSE core functions, DSE applications, and its
practical implementation to support the operation and planning of future power
systems.

  * •

Data Infrastructure for DSE: DSE, whether being implemented in a centralized
or distributed/decentralized manner, needs real-time measurements, and these
measurements need to be transferred to proper locations, such as a control
center, a generation facility, or a substation. The power industry has been
moving forward in the deployment of both sensors and communication networks.
Research topics include: what are the data requirements such as data rates,
signals, sensor placement? What are the data communication requirements, such
as bandwidth, reliability, and redundancy for wide-area system applications?
How to achieve more efficient and robust parallel and/or distributed
implementations of DSE? And for a hierarchical system, how would the data
communication network be structured and where data should be sent to for the
DSE application?

  * •

DSE Core Functions: With the active efforts by many researchers for more than
a decade, many DSE algorithms have been developed, and the DSE performance has
been significantly improved, as summarized in this paper. However, power
systems are evolving to be more complex and different, so DSE core functions
should continue to improve to meet the operational requirements of future
power systems. Research topics include: How to improve computational aspects
of DSE to meet real-time requirements for large-scale systems? How to
formulate DSE with fast inverter dynamics and potentially deal with mixed slow
and fast dynamics? Can machine learning or other data analytical approaches
help improve DSE performance in the context of a large amount of data and a
large number of system configurations? And how to break DSE into pieces in
real-time for islanding situations, especially for a more resilient power
system?

  * •

Development of DSE Applications: DSE provides unprecedented detailed dynamic
information, compared with SSE. A large number of applications in modeling,
monitoring and operation can benefit from DSE as mentioned in this paper. Many
more are emerging. There are conventional applications that can be enhanced by
DSE-provided dynamic information and new applications that can be enabled by
DSE. Some of them include: aggregated model calibration for wind/solar farms,
loads and DERs; oscillation source location; look-ahead DSA; visibility and
detection of converter-induced instability; anomaly detection of DERs, where
the complicated control loop may have fault or failures etc.

  * •

Practicality of DSE Applications: All the research efforts have a common goal
which is to make DSE a practical part of the power system functions. Model
calibration is an example of DSE applications, which is already in commercial
tools, reliability standards, and industry practices. For most DSE
applications, there are still significant gaps to address in terms of their
practicality. Research questions include: How to transition from SSE to DSE
while SSE and DSE will most likely co-exist for a significant period? How to
make DSE compatibility with the control room environment in terms of its
information technology infrastructure and human-machine interface? For
example, DSE results are updated too fast for operators to take actions and
how to present the critical information only to operators would be an
important research need. This may be addressed by development of advanced AI
tools fed with DSE data streaming for visibility and stability assessment;
What training should be developed and provided to prepare the workforce for
DSE applications? And what new standards are needed to enable the DSE
applications?

## References

  * [1] A. Abur, “Dynamic State Estimation,” In Smart Grid Handbook (eds C. Liu, S. McArthur and S. Lee), 2016.
  * [2] B. Kroposki, _et al._ , “Achieving a 100% renewable grid: operating electric power systems with extremely high levels of variable renewable energy,” _IEEE Power and Energy Magazine_ , vol. 15, no. 2, pp. 61-73, 2017.
  * [3] Peter W. Sauer, M. A. Pai, Joe H. Chow. _Power System Dynamics and Stability: With Synchrophasor Measurement and Power System Toolbox_ , 2nd Edition, Wiley, 2019.
  * [4] J. Zhao, _et al._ , “Power system dynamic state estimation: motivations, definitions, methodologies, and future work,” _IEEE Trans. Power Syst._ , vol. 34, no. 4, pp. 3188-3198, July 2019.
  * [5] W. Miller, J. Lewis, “Dynamic state estimation in power systems,” _IEEE Trans. Automatic Control_ , vol. 16, no. 6, pp. 841-846, 1971.
  * [6] J. De La Ree, V. Centeno, J. Thorp, A. Phadke, “Synchronized phasor measurement applications in power systems,” _IEEE Trans. Smart Grid_ , vol. 1, no. 1, pp. 20-27, 2010.
  * [7] E. Song, K. Lee, G. FitzPatrick, Y. Zhang, “Interoperability test for IEC 61850-9-2 standard-based merging units,” In _IEEE Power and Energy Society Innovative Smart Grid Technologies Conference_ , pp. 1-6, 2017.
  * [8] Z. Huang, K. Schneider, J. Nieplocha, “Feasibility studies of applying Kalman filter techniques to power system dynamic state estimation,” in _Proceedings of the 8th International Power Engineering Conference_ , Singapore, 3-6 December 2007.
  * [9] E. Ghahremani and I. Kamwa, “Online state estimation of a synchronous generator using unscented Kalman filter from phasor measurements units,” _IEEE Trans. Energy Convers._ , vol. 26, no. 4, pp. 1099–1108, Dec. 2011.
  * [10] N Zhou, D. Meng, Z. Huang, G. Welch, “Dynamic state estimation of a synchronous machine using PMU data: A comparative study,” _IEEE Trans. Smart Grid_ , vol. 6, no. 1, pp. 450-460, Jan. 2015.
  * [11] M. Netto and L. Mili, “A robust data-driven Koopman Kalman filter for power systems dynamic state estimation,” _IEEE Trans. Power Syst._ , vol. 33, no. 6, pp. 7228-7237, Nov. 2018.
  * [12] A. Rouhani and A. Abur, “Observability analysis for dynamic state estimation of synchronous machines,” _IEEE Trans. Power Syst._ , vol. 32, no. 4, pp. 3168-3175, July 2017.
  * [13] J. Zhao, M. Netto and L. Mili, “A robust iterated extended Kalman filter for power system dynamic state estimation,” _IEEE Trans. Power Syst._ , vol. 32, no. 4, pp. 3205-3216, July 2017.
  * [14] J. Zhao and L. Mili, “A theoretical framework of robust H-infinity unscented Kalman filter and its application to power system dynamic state estimation,” _IEEE Trans. Signal Proc._ , vol. 67, no. 10, pp. 2734-2746, 2019.
  * [15] A. Rouhani and A. Abur, “Linear phasor estimator assisted dynamic state estimation,” _IEEE Trans. Smart Grid_ , vol. 9, no. 1, pp. 211-219, Jan. 2018.
  * [16] L. Wang, K. Morison, “Implementation of online security assessment,” _IEEE Power and Energy Magazine_ , vol. 4, no. 5, pp. 46-59, 2006.
  * [17] Z. Huang, S. Jin, Y. Chen, R. Diao, “Predictive dynamic security assessment through advanced computing,” in _Proceedings of the IEEE Power and Energy Society General Meeting_ , Washington, DC, USA, July 27-31, 2014.
  * [18] IEEE Recommended Practice for Excitation System Models for Power System Stability Studies, IEEE Std 421.5-2016 (Revision of IEEE Std 421.5-2005), vol., no., pp.1-207, 26 Aug. 2016.
  * [19] H. Su, R. Mutukutti, D. Apps, “Impacts of variable quadrature reactance on power system stabilizer performance,” _IEEE Power Engineering Society General Meeting_ , Vancouver, pp. 1-5, 2013.
  * [20] A. K. Singh, B. C. Pal, “Rate of change of frequency estimation for power systems using interpolated DFT and Kalman filter,” _IEEE Trans. Power Syst._ , vol. 34, no. 4, pp. 2509-2517, 2019.
  * [21] J. Qi, K. Sun, W. Kang, “Optimal PMU placement for power system dynamic state estimation by using empirical observability gramian,” _IEEE Trans. Power Syst._ , vol. 30, no. 4, pp. 2041–2054, Jul. 2015.
  * [22] A. Abur and A. Gomez-Exposito, _Power system state estimation: theory and application_. New York, USA: Marcel Dekker, 2004.
  * [23] N. Zhou, S. Wang, J. B. Zhao, Z. Huang, “Application of detectability analysis for power system dynamic state estimation”, _IEEE Trans. Power Syst._ , 2020.
  * [24] M. Gol and A. Abur, “LAV based robust state estimation for systems measured by PMUs,” _IEEE Trans. Smart Grid_ , vol. 5, no. 4, pp. 1808– 1814, Apr. 2014.
  * [25] A. Rouhani and A. Abur, “A robust dynamic state estimator against exciter failures,” in _Proc. North Amer. Power Symp._ , pp. 1–6, 2016.
  * [26] G. Anagnostou, _et. al_ , “Observer-based anomaly detection of synchronous generators for power systems monitoring,” _IEEE Trans. Power Syst._ , vol. 33, no. 4, pp. 4228–4237, Jul. 2018.
  * [27] S. Wei, _et. al_ , “Model-free MLE estimation for online rotor angle stability assessment with PMU data,” _IEEE Trans. Power Syst._ , vol. 33, no. 3, pp. 2463-2476, 2018.
  * [28] A. P. S. Meliopoulos, _et al._ , “Dynamic state estimation-based protection: status and promise,” _IEEE Trans. Power Delivery_ , vol. 32, no. 1, pp. 320-330, Feb. 2017.
  * [29] J. B. Zhao, L. Mili, “A robust generalized-maximum likelihood unscented Kalman filter for power system dynamic state estimation,” _IEEE Journal of Selected Topics in Signal Processing_ , vol. 12, no. 4, pp. 578-592, 2018.
  * [30] R. D. Quint, “Power plant model verification guideline,” 2016.
  * [31] Z. Huang, P. Du, D. Kosterev, S. Yang, “Generator dynamic model validation and parameter calibration using phasor measurements at the point of connection,” _IEEE Trans. Power Syst._ , vol. 28, no. 2, pp.1939-1949, May 2013.
  * [32] R. Huang, _et al._ , “Calibrating parameters of power system stability models using advanced ensemble Kalman filter,” _IEEE Trans. Power Syst._ , vol. 33, no. 3, pp. 2895-2905, May 2018.
  * [33] S. Akhlaghi, S. Raheem, N. Zhou, “Model validation lessons learned through implementing NERC MOD-033-1,” in _Proc. of IEEE PES General Meeting_ , Montreal, CA, Aug. 2020.
  * [34] Y. Zhang, E. Muljadi, D. Kosterev, M. Singh, “Wind power plant model validation using synchrophasor measurements at the point of interconnection,” _IEEE Trans. Sustainable Energy_ , vol. 6, no. 3, pp. 984-992, 2015.
  * [35] A. Rouhani and A. Abur, “Constrained iterated unscented Kalman filter for dynamic state and parameter estimation,” _IEEE Trans. Power Syst._ , vol. 33, no. 3, pp. 2404–2414, May 2018.
  * [36] G. Valverde, E. Kyriakides, G. Heydt, V. Terzija “Nonlinear estimation of synchronous machine parameters using operating data,” _IEEE Trans. Energy Conversion_ , vol. 26, no. 3, pp. 831-839, 2011.
  * [37] R. Fan, _et. al_ , “Dynamic state estimation and parameter calibration of a DFIG using the ensemble Kalman filter,” in _Proc. IEEE PES General Meeting_ , Denver, CO, USA, Jul. pp. 1–5, 2015.
  * [38] M. Gonzalez-Cagigal, J. Rosendo-Macias, A. Gomez-Exposito, “Parameter estimation of fully regulated synchronous generators using unscented Kalman filters,” _Electric Power Systems Research_ , vol. 168, pp. 210-217, 2019.
  * [39] S. Joo, C. Liu, L. Jones, J. Choe, “A coherency and aggregation techniques incorporating rotor and voltage dynamics,” _IEEE Trans. Power Syst._ , vol. 19. No. 2, pp. 1068-1075, 2004.
  * [40] M. Aghamohammadi, S. Tabandeh, ”A new approach for online coherency identification in power systems based on correlation characteristics of generators rotor oscillations,” _International Journal of Electrical Power & Energy Systems_, vol. 83, pp. 470-484, 2016.
  * [41] A. Chakrabortty, J. H Chow, A. Salazar, “A measurement-based framework for dynamic equivalencing of large power systems using wide-area phasor measurements,” _IEEE Trans. Smart Grid_ , vol. 2, no. 1, 68-81, 2011.
  * [42] A. Rouhani and A. Abur, “Real-time dynamic parameter es-timation for an exponential dynamic load model,” _IEEE Trans. Smart Grid_ , vol. 7, no. 3, pp. 1530–1536, May 2016.
  * [43] M. Arif, B. C. Pal, A. K. Singh, “Estimating dynamic model parameters for adaptive protection and control in power system,” _IEEE Trans. Power Syst._ , vol. 30, no. 2, pp. 829–839, Mar. 2015.
  * [44] P. Bhui, N. Senroy, A. K. Singh, B. C. Pal, “Estimation of inherent governor dead-band and regulation using unscented Kalman filter,” _IEEE Trans. Power Syst._ , vol. 33, no. 4, pp. 3546-3558, 2018.
  * [45] J. Sanchez-Gasca, _et al._ , “Identification of electromechanical modes in power systems,” IEEE Task Force Rep. no. IEEE PES-TR15, Jun. 2012.
  * [46] S. Maslennikov, B. Wang, Q. Zhang, E. Litvinov, “A test cases library for methods locating the sources of sustained oscillations,” In _Proc. IEEE PES General Meeting_ , pp. 1-5, 2016.
  * [47] J. B. Zhao, L. Mili, F. Milano, “Robust frequency divider for power system online monitoring and controls,” _IEEE Trans. Power Syst._ , vol. 33, no. 4, pp. 4414–4423, Jul. 2018.
  * [48] J. B. Zhao, Y. Tang, V. Terzija, “Robust online estimation of power system center of inertia frequency,” _IEEE Trans. Power Syst._ , vol. 34, no. 1, pp. 821–825, Jan. 2019.
  * [49] F. Milano and A. Ortega, “Frequency divider,” _IEEE Trans. Power Syst._ , vol. 32, no. 2, pp. 1493–1501, Mar. 2017.
  * [50] S. Wang, J. Zhao, Z. Huang, R. Diao, “Assessing Gaussian assumption of PMU measurement error using field data,” _IEEE Trans. Power Del._ , vol. 33, no. 6, pp. 3233–3236, Dec. 2017.
  * [51] J. B. Zhao, L. Mili, A. Gomez-Exposito, “Constrained robust unscented Kalman filter for generalized dynamic state estimation,” _IEEE Trans. Power Syst._ , vol. 34, no. 5, pp. 3637-3646, 2019.
  * [52] Y. Chakhchoukh, H. Lei, B. Johnson, “Diagnosis of outliers and cyber attacks in dynamic PMU-based power state estimation,” _IEEE Trans. Power Syst._ , vol. 35, no. 2, pp. 1188-1197, 2020.
  * [53] Q. Yang, T. Bi, J. Wu, “WAMS implementation in China and the challenges for bulk power system protection,” in _Proc. IEEE Power Eng. Soc. General Meeting_ , Jun, pp. 1–6, 2007.
  * [54] G. Anagnostou, B. Pal, “Impact of overexcitation limiters on the power system stability margin under stressed conditions,” _IEEE Trans. Power Syst._ , vol. 31, no. 3, pp. 2327–2337, May 2016.
  * [55] E. Ghahremani, I. Kamwa, “Dynamic state estimation in power system by applying the extended Kalman filter with unknown inputs to phasor measurements,” _IEEE Trans. Power Syst._ , vol. 26, no. 4, pp. 2556-2566, Nov. 2011.
  * [56] E. Ghahremani, I. Kamwa, “Local and wide-area PMU-based decentralized dynamic state estimation in multi-machine power systems,” _IEEE Trans. Power Syst._ , vol. 31, no. 1, pp. 547-562, 2016.
  * [57] G. Anagnostou, B. Pal, “Derivative-free Kalman filtering based approaches to dynamic state estimation for power systems with unknown inputs,” _IEEE Trans. Power Syst._ , vol. 33, no. 1, pp. 116-130, 2017.
  * [58] J. B. Zhao, _et. al_ , “Correlation-aided robust decentralized dynamic state estimation of power systems with unknown control inputs,” _IEEE Trans. Power Syst._ , 2019.
  * [59] P. Marchi, P. Estevez, F. Messina, C. Galarza, “Loss of excitation detection in synchronous generators based on dynamic state estimation,” _IEEE Trans. Energy Conversion_ , 2020.
  * [60] M. Diaz-Aguilo and F. de Leon, “Adaptive soil model for real-time thermal rating of underground power cables,” _IET Science, Measurement & Technology_, vol. 9, no. 6, pp. 654-660, 9 2015.
  * [61] A. Gomez-Exposito, _et. al_ , “City-friendly smart network technologies and infrastructures: the spanish experience,”  _Proceedings of the IEEE_ , vol. 106, no. 4, pp. 626-660, 2018.
  * [62] I. Ortega Rivera, _et. al_ , “A dynamic state estimator based control for power system damping,” _IEEE Trans. Power Syst._ , vol. 33, no. 6, pp. 6839-6848, 2018.
  * [63] A. Paul, I. Kamwa, G. Joos, “PMU signals responses-based RAS for instability mitigation through on-the fly identification and shedding of the run-away generators,” _IEEE Trans. Power Syst._ , 2019.
  * [64] “Stability definitions and characterization of dynamic behavior in systems with high penetration of power electronic interfaced technologies,” IEEE Technical Report, IEEE Power and Energy Society, Power System Dynamic Performance Committee, 2020.
  * [65] S. Swain, B. Subudhi, “Grid synchronization of a PV system with power quality disturbances using unscented Kalman filtering,” _IEEE Trans. Sustainable Energy_ , vol. 10, no. 3, pp. 1240-1247, 2018.
  * [66] S. Yu, _et. al_ , “State estimation of doubly fed induction generator wind turbine in complex power systems,” _IEEE Trans. Power Syst._ , vol. 31, no. 6, pp. 4935-4944, 2016.
  * [67] S. Yu, T. Fernando, H. H.-C. Iu, K. Emami, “Realization of state-estimation-based DFIG wind turbine control design in hybrid power systems using stochastic filtering approaches,” _IEEE Trans. Industrial Informatics_ , vol. 12, no. 3, pp. 1084-1092, 2016.
  * [68] S. Yu, _et. al_ , “An unscented particle filtering approach to decentralized dynamic state estimation for DFIG wind turbines in multi-area power systems,” _IEEE Trans. Power Syst._ , 2019.
  * [69] G. Anagnostou, L. Kunjumuhammed, B. Pal, “Dynamic state estimation for wind turbine models with unknown wind velocity,” _IEEE Trans. Power Syst._ , vol. 34, no. 5, pp. 3879-3890, 2019.
  * [70] M. A. Gonzalez-Cagigal, J. A. Rosendo-Macias, A. Gomez-Exposito, “Parameter estimation of wind turbines with PMSM using cubature Kalman filters,” _IEEE Trans. Power Syst._ , vol. 35, no. 3, pp. 1796-1804, May 2020.
  * [71] “1,200 MW Fault Induced Solar Photovoltaic Resource Interruption Disturbance Report: Southern California 8/16/2016 Event,” NERC, Atlanta, GA, June 2017.


††thanks: Corresponding author<EMAIL_ADDRESS>Corresponding author<EMAIL_ADDRESS>
# Autopoietic influence hierarchies in pancreatic $\beta$ cells

Dean Korošak  Institute of Physiology, Faculty of Medicine, University of
Maribor, 2000 Maribor, Slovenia  Faculty of Civil Engineering, Transportation
Engineering and Architecture, University of Maribor, 2000 Maribor, Slovenia
Marko Jusup  Tokyo Tech World Research Hub Initiative (WRHI), Institute of
Innovative Research, Tokyo Institute of Technology, Tokyo 152-8552, Japan
Boris Podobnik  Faculty of Civil Engineering, University of Rijeka, 51000
Rijeka, Croatia  Center for Polymer Studies, Boston University, Boston MA
02215, USA  Zagreb School of Economics and Management, 10000 Zagreb, Croatia
Luxembourg School of Business, 2453 Luxembourg, Luxembourg  Faculty of
Information Studies in Novo mesto, 8000 Novo Mesto, Slovenia  Andraž Stožer
Institute of Physiology, Faculty of Medicine, University of Maribor, 2000
Maribor, Slovenia  Jurij Dolenšek  Institute of Physiology, Faculty of
Medicine, University of Maribor, 2000 Maribor, Slovenia  Faculty of Natural
Sciences and Mathematics, University of Maribor, 2000 Maribor, Slovenia
Petter Holme  Tokyo Tech World Research Hub Initiative (WRHI), Institute of
Innovative Research, Tokyo Institute of Technology, Tokyo 152-8552, Japan
Marjan Slak Rupnik  Institute of Physiology, Faculty of Medicine, University
of Maribor, 2000 Maribor, Slovenia  Center for Physiology and Pharmacology,
Medical University of Vienna, 1090 Vienna, Austria  Alma Mater Europaea –
European Center Maribor, 2000 Maribor, Slovenia

(August 27, 2024)

###### Abstract

$\beta$ cells are biologically essential for humans and other vertebrates.
Because their functionality arises from cell-cell interactions, they are also
a model system for collective organization among cells. There are currently
two contradictory pictures of this organization: the hub-cell idea pointing at
leaders who coordinate the others, and the electrophysiological theory
describing all cells as equal. We use new data and computational modeling to
reconcile these pictures. We find via a network representation of interacting
$\beta$ cells that leaders emerge naturally (confirming the hub-cell idea),
yet all cells can take the hub role following a perturbation (in line with
electrophysiology).

The importance of the nutrient-sensing and insulin-secreting $\beta$ cells in
vertebrates is hard to overstate. These cells reside in pancreatic islets,
where they extensively communicate with each other and their environment [1].
The intercellular communication, in particular, serves to coordinate and
synchronize cellular operations through which insulin is released in
proportion to stimulation and metabolic requirements [2, 3]. The delicate
nature of this task is reflected in the continual need to prevent
oversecretion, and subsequent hypoglycemia [4], despite intracellular stores
holding sufficient insulin to exceed the lethal dose by orders of magnitude if
released at once.

On a microscopic scale, $\beta$ cells are electrically and metabolically
coupled via gap junctions built from the protein connexin-36 (Cx36) [5]. The
key role of such coupling is seen in the fact that the loss of Cx36 channels
has a detrimental impact on $\beta$-cell cooperation [6, 7], leading to
uncoordinated plasma depolarizations, the desynchronization of calcium
signals, and increases in basal insulin release. Sufficient Cx36 coupling, by
contrast, curbs the $\beta$-cell intrinsic heterogeneity in glucose
sensitivity [8], and in mouse islets, limits the threshold for stimulatory
glucose concentration [7] to a narrow band around $\approx$7 mM.
Metaphorically, coupled $\beta$ cells are like individual soldiers who fall in
line when the communication channels between them are open. A typical $\beta$
cell is coupled to between six and eight neighbors [9, 10]. According to
electrophysiological theory [11], single cells lack mechanisms to become
pacemakers beyond their immediate neighborhood. In terms of our metaphor,
there are no generals in the army.

The above-mentioned lateral organization is diametrically opposed to the
picture by functional studies of a large number of communicating $\beta$ cells
(referred to as $\beta$-cell collectives). For example, intrinsic cellular
heterogeneity in glucose sensing [8], heterogeneous gap-junction coupling
[12], and extensive paracrine signaling [13] all contribute to an islet-wide
complicated cytoplasmic Ca2+ dynamics. This dynamics of the Ca2+ concentration
in the cytosol [14, 15] is a key insulin-secretion driver [16]. It has
recently become possible to record the Ca2+ dynamics with a great spatial and
temporal precision using the functional multicellular confocal imaging [17,
18]. The high data content of such imaging enables mapping the functional
organization of a $\beta$-cell collective onto a complex network such that a
pair of cells (i.e., network nodes) is linked if the cross-correlation between
the corresponding Ca2+ signals is large enough [19, 20, 21, 22, 23, 24].

Complex-network representations of the functional organization of $\beta$-cell
collectives (i.e., functional networks) that arise from calcium cross-
correlations in intact pancreatic islets reveal a modular structure [25]
intertwined with small-world [26, 20] properties. Of particular interest is
that some studies [22, 24] point to small subsets of highly active $\beta$
cells whose connectedness and impact on synchronization across islets make
these cells candidate leaders. Because of their large degree, we call such
candidate leaders hub cells. The proponents of the hub-cell idea list many
reasons [27] why electrophysiological measurements may have missed detecting
hubs. Overall, functional networks suggest that our army of $\beta$ cells is
led by a few generals.

We will hereafter show that the hub-cell idea, rather than contradicting,
actually complements electrophysiology. Based on a remarkable agreement
between empirical and computational findings, we will argue in favor of an
autopoietic influence hierarchy among $\beta$ cells. Leader cells emerge
through cell-cell communication, and thus need not be genetically predisposed
for leadership. In the language of the army metaphor, generals do lead, but
not by birthright; they get promoted by their peers. Leader cells,
furthermore, remain under the radar of electrophysiological measurements by
communicating with immediate neighbors just as any other cell would. A direct
implication is that a more-or-less arbitrary cell could emerge as a leader,
which in turn makes for an ultra robust architecture of $\beta$-cell
collectives. Additional wide-reaching implications underpin some of the vital
$\beta$-cell features.

##  I Methods

We based our empirical analyses on a dataset obtained via Ca2+ imaging of an
acutely prepared pancreatic tissue slice [28, 17] comprising a rodent oval-
shaped islet (approx. dimensions: 370 µm$\times$200 µm). Ca2+ signals were
recorded using a functional multi-cellular imaging technique at 10 Hz and
256$\times$256 pixel resolution in 8-bit grayscale color depth. The freely
downloadable dataset [29] consisted of 65,536 Ca2+ signals, each with 23,873
data points. Our focus was on fast oscillations, which is why all signals were
detrended and standardized before use [30].

For the purpose of constructing empirical functional networks, we randomly
picked $N=100$ signals, denoted $x_{i}(t)$, from the dataset and computed
cross-correlations $c_{ij}=[x_{i}x_{j}]$, where $[\cdot]$ is a time-averaging
operator. The number of data points in each signal was $L=3\times 10^{3}$,
corresponding to 5-minute exposures of the pancreatic tissue slice to the
stimulating glucose concentration of 8 mM. We link two signals $i$ and $j$ in
the network if $c_{ij}>c_{0}$, where $c_{0}$ is a threshold selected to give a
mean degree $\langle{}k\rangle=10$. We have previously shown [31] that the
bulk cross-correlation spectrum of typical Ca2+ signals follows the
predictions of random matrix theory, but the states outside the bulk-spectrum
edges carry biological information. The delocalized states corresponding to
the largest eigenvalues of the cross-correlation matrix, in particular, were
found to harbor contributions from all Ca2+ signals, thus revealing an islet-
wide collective mode [32].

Figure 1:  Functional networks of $\beta$ cells contain high-degree hubs. The
plot shows the configuration-ensemble degree distribution, $P$, as a function
of the node degree, $k$.

The computational aspects of the study comprised constructing functional
networks from simulated Ca2+ signals. Simulations were based on our model [30]
that had been found to mimic empirical signals closely. Compared to typical
electrophysiological modeling by means of differential-algebraic systems [21],
the model’s structure is very simple. Each of the $N=100$ nodes represents a
$\beta$ cell arranged in a random regular network with the degree distribution
$P(k)=\delta(k-k_{0})$ such that trivially the average node degree is
$\langle{}k\rangle=k_{0}$. Nodes can change their binary state, indicating the
calcium activity of individual $\beta$ cells, from active to inactive or vice
versa in two ways. Internal activation is controlled by a forcing parameter,
$f$, that is interpretable as the inverse of the glucose concentration to
which $\beta$ cells are exposed. External activation is controlled by the
state of the $k_{0}$ nearest neighbors in conjunction with the coupling
strength, $r$, that is interpretable as the intensity of gap-junctional ion
exchange. Functional networks were constructed from simulated Ca2+ signals in
exactly the same way as from empirical ones.

##  II Results

The degree distributions of empirical and simulated functional networks are
practically the same (Fig. 1). We distinguish between two cases. Open symbols
pertain to functional networks extracted from empirical data per se (squares),
and from simulations with strong coupling, $h=0.8$ (circles). In this case,
the obtained degree distributions closely follow a negative binomial
distribution $P(k)=\binom{k+r-1}{k}p^{k}(1-p)^{r}$, where
$p=\langle{}k\rangle/(r+\langle{}k\rangle)\approx{}0.77$ and $r=3$ is a
numerical parameter. The opposite case, shown using filled symbols, consists
of functional networks extracted upon randomizing empirical data (squares) and
running simulations with weak coupling, $h=0.1$ (circles). Both data
randomization and weak coupling should produce random functional networks,
what indeed transpires given that the obtained degree distributions closely
follow a binomial distribution $P(k)=\binom{N}{k}p^{k}(1-p)^{N-k}$ with
$p=\langle{}k\rangle/N$. The binomial distribution in turn converges to a
Poisson distribution $P(k)=\langle{}k\rangle^{k}\exp(-\langle{}k\rangle)/k!$
in the large network limit $N\to\infty$ when $\langle{}k\rangle$ is fixed.

Figure 2:  $\beta$-cell functional networks exhibit positive assortative
mixing. The plot compares the elements $a_{ij}$ of the configuration-ensemble
adjacency matrix against the normalized degree product,
$(k_{i}k_{j}/\langle{}k\rangle N)^{1+\beta}$, for two different values of the
$\beta$ exponent, $\beta=0$ and $\beta=0.1$.

Aside from degree distributions, another informative way of characterizing
functional networks is degree correlations, also known as degree
assortativity. Ref. [33] explains that the elements $a_{ij}$ of a
configuration-ensemble adjacency matrix should satisfy the strictly linear
relationship $a_{ij}=k_{i}k_{j}/\langle{}k\rangle N$ in zero-assortativity
situations. Otherwise, $a_{ij}=(k_{i}k_{j}/\langle{}k\rangle N)^{1+\beta}$,
where $\beta<0$ ($\beta>0$) indicates negative (positive) assortative mixing.
Plots of $a_{ij}$ against $k_{i}k_{j}/\langle{}k\rangle N$ (Fig. 2), while
once again revealing an agreement between empirical and simulated functional
networks, primarily show that $\beta\approx 0.1$. The degree assortativity of
the functional networks is therefore positive, meaning that network nodes are
preferentially linked to other nodes of similar degree.

Figure 3:  Assortativity and the largest adjacency-matrix eigenvalue of
functional networks are robust to random perturbations, but not to
perturbations targeting hubs. The plot shows assortativity, $\rho$, and the
largest eigenvalue of the configuration-ensemble adjacency matrix,
$\lambda_{1}$, (from which the average node degree, $\langle{}k\rangle$, has
been subtracted) both as the functions of the fraction of perturbed nodes,
$N_{\mathrm{r}}/N$. Figure 4:  Hub cells carry most of the cross-correlation
content of a $\beta$-cell collective. (a) Empirical neighborhood-wide average
cross-correlation of nodes, $c$, as a function of the node degree, $k$. Open
circles, squares, and triangles respectively denote no perturbation, random
perturbation, and hub randomization. (b) Same as (a), but simulated. Open
circles and triangles respectively denote no perturbation and random
perturbation when coupling, $h=0.8$ is strong, whereas filled circles denote
no perturbation when coupling, $h=0.1$, is weak. (c) The simulated average
number of active neighbors, $k^{*}$, as a function of the node degree. Open
and filled circles respectively denote the strong, $h=0.8$, and the weak
coupling, $h=0.1$.

The results so far allow us to estimate the assortativity coefficient of
functional networks in an alternative way, which we can then compare with the
definition [34, 35], that is, with the degree correlation coefficient between
pairs of linked nodes
$\rho=(\langle{}k\rangle\langle{}k^{2}k_{\mathrm{nn}}\rangle-\langle{}k^{2}\rangle^{2})/(\langle{}k\rangle\langle{}k^{3}\rangle-\langle{}k^{2}\rangle^{2})$,
where $k_{\mathrm{nn}}^{i}=k_{i}^{-1}\sum_{j}a_{ij}k_{j}$ is the average
degree of the $i$th node’s nearest neighbors. Using the exponent $\beta\approx
0.1\ll 1$ and the fact that the aforementioned negative binomial distribution
converges to a gamma distribution $P(k)=q^{-r}k^{r-1}\exp(-k/q)/\Gamma(r)$ for
large $q=\langle{}k\rangle/r$, the assortativity coefficient is
$\rho\approx\beta{}q^{2\beta}\approx{}0.11$ (see Ref. [33]). Following the
definition of $\rho$ indeed yields a very similar value (Fig. 3; open symbols
for $N_{\mathrm{r}}/N=0$).

We furthermore tested how two different perturbations [36], signal
randomization and replacement, affect the assortativity coefficient. Here,
randomization and replacement respectively mean that a fraction
$N_{\mathrm{r}}/N$ of signals used in the construction of the original
configuration ensemble were either randomly reshuffled or replaced with other
random signals from the broader dataset. We found that assortativity of
functional networks is robust to random perturbations, but not to
perturbations targeting hubs (Fig. 3; open symbols for $N_{\mathrm{r}}/N>0$).
It is therefore hubs that preferentially link to other hubs, which emphasizes
their importance in the overall structure of functional networks.
Specifically, we are seeing the evidence of a purposeful correlating
mechanisms at work among $\beta$-cells because otherwise heterogeneous
networks are expected to be disassortative [33, 35].

Another measure characterizing a network’s structure is the largest eigenvalue
of the adjacency matrix, $\lambda_{1}$. Similarly as with the assortativity
coefficient, the value of $\lambda_{1}$ is robust to random perturbations, but
not to perturbations targeting hubs (Fig. 3). As progressively more hubs get
perturbed, the maximum adjacency-matrix eigenvalue falls to
$\lambda_{1}=1+\langle{}k\rangle$, which is the value characteristic of random
networks with a Poisson degree distribution. This last result reaffirms the
importance of hubs for the structure of functional networks and, in turn,
hints at an interesting connection with recent experimental findings that
generated much excitement.

Experiments show [22, 24] that perturbing hub cells can quickly decorrelate a
$\beta$-cell collective. This implies that hubs carry most of the collective’s
cross-correlation content, which could be tested using our empirical and
simulated Ca2+ signals alike. To this end, we employed the pairwise cross-
correlations $c_{ij}$ to define $c_{i}=k_{i}^{-1}{\sum_{j}a_{ij}c_{ij}}$, that
is, the $i$th node’s neighborhood-wide average cross-correlation. By pairing
the $c_{i}$ values with the corresponding node degree $k_{i}$, we obtained a
function $c(k)$. If the $\beta$-cell collective’s cross-correlation content
were equally distributed among all nodes, then we should see
$c(k)=\mathrm{const.}$, whereas if hubs carried most of the content, then
$c(k)$ should be an increasing function of $k$. We find that the latter is
true both for empirical (Fig. 4a) and simulated (Fig. 4b) Ca2+ signals. As
with assortativity $\rho$ and the largest eigenvalue $\lambda_{1}$, random
perturbations have little bearing on the function $c(k)$, but targeting high-
degree nodes causes decorrelation. With the help of the dynamical network
model, we see that perturbing about 10 % of hubs is very similar to running
simulations with a weak coupling of $h=0.1$ (Fig. 4b).

The model additionally allows us to examine whether the node-degree dependence
of the function $c(k)$ could be due to hubs somehow communicating with more of
their neighbors than lower-degree nodes. The number of active neighbors
$k^{*}$ is, as expected, independent of the degree $k$ when the coupling is
weak, $h=0.1$ (Fig. 4c). The same, however, approximately holds even when the
coupling is strong, $h=0.8$ (Fig. 4c). Because, in the model, all nodes are
equal and none of them communicate with an unusual number of neighbors, we are
left with a conclusion that the emergence of hubs in functional networks is
endogenous to cell-cell communication.

##  III Discussion

Herein, we used a combination of empirical and computational methods to shed
new light on, among others, a fundamental tension between electrophysiological
theory and the hub-cell idea as it pertains to pancreatic $\beta$ cells. We
started by constructing functional $\beta$-cell networks from empirical Ca2+
signals, and after that estimated the configuration ensemble [37, 38] of such
networks. Meanwhile, we also simulated fast $\beta$-cell activity using a
dynamical network model [39, 40] whose outputs faithfully mimic the properties
of said empirical Ca2+ signals [30]. Upon repeating the construction of
functional networks, but now from simulated signals, we found remarkable
quantitative agreement between the results.

Irrespective of the signal origin, empirical or simulated, the properties of
functional networks—such as assortativity or the largest eigenvalue of the
configuration-ensemble adjacency matrix—all critically depend on the presence
of highly connected (i.e., hub) nodes. Our results thus support the hub-cell
idea. A similar heterogeneous organization of functional networks has also
been discovered in the collectives of chemosensing cells [41, 42]
communicating via gap junctions.

The support for the hub-cell idea would be a strong blow to
electrophysiological theory if nodes in the dynamical network model had any
distinctive properties, for example, different degrees or an internal
structure. In the model, however, all nodes are completely identical,
upholding the ideas of electrophysiological theory that no cell is predisposed
for leadership. Leader cells, in fact, emerge through cell-cell communication
as evidenced by tracing the degree distributions of functional networks to the
coupling strength in the dynamical network model. In the language of the army
metaphor, generals do lead, but not by birthright; they get promoted by their
peers. Our findings thus not only reconcile the concept of hubs with
electrophysiological theory, but also point to an extremely robust
architecture of collectives.

Before expanding on the last remark about the robustness of $\beta$-cell
collectives, let us examine a physiological advantage provided by the
heterogeneity of functional networks. Cells have been shown to use cell-cell
distance to optimize their sensing precision [43]. In compact microorgans such
as pancreatic islets, cells cannot easily adjust mutual distances, but could
instead use functional-network heterogeneity to sharpen their collective
response. Let $\lambda_{1}$ be the largest adjacency-matrix eigenvalue and
$\mathbf{v}^{1}$ the corresponding eigenvector [44, 32]. Because the
components of the latter vector quantify the importance of network nodes, we
assume that
$y_{i}=\sum_{j\in\Omega_{i}}v_{j}^{1}x_{j}=\sum_{j}a_{ij}v_{j}^{1}x_{j}$
represents the signal that node $i$ integrates from the nearest neighbors
$j\in\Omega_{i}$, while all other signals are treated as constant-variance
noise. Under this assumption, the signal-to-noise ratio of sensing cells is
proportional to the variance
$\operatorname{\mathbb{E}}[y^{2}]=\sum_{ijk}a_{ij}v_{j}^{1}a_{ik}v_{k}^{1}\operatorname{Cov}(x_{j},x_{k})\gtrapprox
c_{0}\lambda_{1}^{2}$, where $c_{0}$ is the cross-correlation threshold used
during the construction of functional networks. The sensing precision is thus
proportional to $\lambda_{1}^{2}$, which would equal
$\lambda_{1}=1+\langle{}k\rangle$ if functional networks were random, but
increases to $\lambda_{1}=1+\langle{}k\rangle+\langle{}k\rangle/r$ because of
heterogeneity. Accordingly, cell-cell communication imparts a sharp glucose-
sensing acuity to $\beta$-cell collectives.

The example on glucose-sensing acuity demonstrates a physiological advantage
of heterogeneous functional networks over homogeneous ones. This advantage
would, however, come at a cost of vulnerability to hub-node failures [36] if
heterogeneity were imprinted in the underlying physical network of $\beta$
cells. Our results instead strongly favor an interpretation by which
influential cells materialize endogenously within $\beta$-cell collectives,
giving rise to influence hierarchies that are autopoietic, both in the self-
producing and self-maintaining sense of this term. All that is needed for hub
cells to emerge among identical peers is cell-cell communication, and should a
hub fail, there is no obstacle for another one to re-emerge as long as
communication remains feasible. It would thus seem that in $\beta$-cell
collectives, biology has managed to create an ultra-robust architecture that
is physiologically advantageous as well. This is an impressive feat that could
perhaps inspire thinking about the design and engineering of next-generation
critical infrastructure.

– – – – –

Acknowledgements. DK, AS, and MSR received financial support from the
Slovenian Research Agency (research core funding program no. P3–0396 and
projects no. N3-0048, no. N3-0133 and no. J3-9289). MSR was further supported
by the Austrian Science Fund / Fonds zur Förderung der Wissenschaftlichen
Forschung (bilateral grants I3562–B27 and I4319–B30). MJ received financial
supported from the Japan Society for the Promotion of Science (JSPS) KAKENHI
(grant no. JP 20H04288) as a co-investigator. BP received financial support
from the Slovenian Research Agency (project no. J5-8236), the University of
Zagreb’s project Advanced methods and technologies in Data Science and
Cooperative Systems (DATACROSS; ref. KK.01.1.1.01.009), and the University of
Rijeka. PH received financial support from JSPS KAKENHI (grant no. JP
18H01655).

Author contributions. All authors contributed substantially to all aspects of
the study.

Conflict of interest. The authors declare no conflict of interest, financial
or otherwise.

## References

  * Pipeleers [1987] D. Pipeleers, Diabetologia 30, 277 (1987).
  * Pipeleers  _et al._ [1982] D. Pipeleers, E. Maes, M. Van De Winkel, _et al._ , Proc. Natl. Acad. Sci. USA 79, 7322 (1982).
  * Korošak and Slak Rupnik [2018] D. Korošak and M. Slak Rupnik, Front. Physiol. 9, 31 (2018).
  * Kolb  _et al._ [2020] H. Kolb, K. Kempf, M. Röhling, and S. Martin, BMC Med. 18, 224 (2020).
  * Meda [2018] P. Meda, Biochim. Biophys. Acta Biomembr. 1860, 124 (2018).
  * Ravier  _et al._ [2005] M. A. Ravier, M. Güldenagel, A. Charollais, A. Gjinovci, D. Caille, G. Söhl, C. B. Wollheim, K. Willecke, J.-C. Henquin, and P. Meda, Diabetes 54, 1798 (2005).
  * Speier  _et al._ [2007] S. Speier, A. Gjinovci, A. Charollais, P. Meda, and M. Rupnik, Diabetes 56, 1078 (2007).
  * Benninger and Hodson [2018] R. K. Benninger and D. J. Hodson, Diabetes 67, 537 (2018).
  * Zhang  _et al._ [2008] Q. Zhang, J. Galvanovskis, F. Abdulkader, C. J. Partridge, S. O. Göpel, L. Eliasson, and P. Rorsman, Philos. Trans. R. Soc. A 366, 3503 (2008).
  * Skelin Klemen  _et al._ [2017] M. Skelin Klemen, J. Dolenšek, M. Slak Rupnik, and A. Stožer, Islets 9, 109 (2017).
  * Satin  _et al._ [2020] L. S. Satin, Q. Zhang, and P. Rorsman, Diabetes 69, 830 (2020).
  * Farnsworth  _et al._ [2014] N. L. Farnsworth, A. Hemmati, M. Pozzoli, and R. K. Benninger, J. Physiol. 592, 4431 (2014).
  * Caicedo [2013] A. Caicedo, in  _Seminars in cell & developmental biology_, Vol. 24(1), edited by S. D. Roper (Elsevier, 2013) pp. 11–21.
  * Berridge  _et al._ [2000] M. J. Berridge, P. Lipp, and M. D. Bootman, Nat. Rev. Mol. Cell Biol. 1, 11 (2000).
  * Colecraft [2020] H. M. Colecraft, Biophys. J. 119, 1472 (2020).
  * Idevall-Hagren and Tengholm [2020] O. Idevall-Hagren and A. Tengholm, in  _Seminars in cell & developmental biology_, Vol. 103, edited by W. Han and D. Eizirik (Elsevier, 2020) pp. 20–30.
  * Stožer  _et al._ [2013a] A. Stožer, J. Dolenšek, and M. S. Rupnik, PLOS One 8, e54638 (2013a).
  * Dolenšek  _et al._ [2020] J. Dolenšek, M. S. Klemen, M. Gosak, L. Križančić-Bombek, V. Pohorec, M. S. Rupnik, and A. Stožer (2020), bioRxiv 2020.03.11.986893.
  * Hodson  _et al._ [2012] D. J. Hodson, M. Schaeffer, N. Romanò, P. Fontanaud, C. Lafont, J. Birkenstock, F. Molino, H. Christian, J. Lockey, D. Carmignac, _et al._ , Nat. Commun. 3, 605 (2012).
  * Stožer  _et al._ [2013b] A. Stožer, M. Gosak, J. Dolenšek, M. Perc, M. Marhl, M. S. Rupnik, and D. Korošak, PLOS Comput. Biol. 9, e1002923 (2013b).
  * Cherubini  _et al._ [2015] C. Cherubini, S. Filippi, A. Gizzi, and A. Loppini, Phys. Rev. E 92, 042702 (2015).
  * Johnston  _et al._ [2016] N. R. Johnston, R. K. Mitchell, E. Haythorne, M. P. Pessoa, F. Semplici, J. Ferrer, L. Piemonti, P. Marchetti, M. Bugliani, D. Bosco, _et al._ , Cell Metab. 24, 389 (2016).
  * Gosak  _et al._ [2018] M. Gosak, R. Markovič, J. Dolenšek, M. S. Rupnik, M. Marhl, A. Stožer, and M. Perc, Phys. Life Rev. 24, 118 (2018).
  * Salem  _et al._ [2019] V. Salem, L. D. Silva, K. Suba, E. Georgiadou, S. N. M. Gharavy, N. Akhtar, A. Martin-Alonso, D. C. Gaboriau, S. M. Rothery, T. Stylianides, _et al._ , Nat. Metab. 1, 615 (2019).
  * Markovič  _et al._ [2015] R. Markovič, A. Stožer, M. Gosak, J. Dolenšek, M. Marhl, and M. S. Rupnik, Sci. Rep. 5, 7845 (2015).
  * Rutter and Hodson [2013] G. A. Rutter and D. J. Hodson, Mol. Endocrinol. 27, 1984 (2013).
  * Rutter  _et al._ [2020] G. A. Rutter, N. Ninov, V. Salem, and D. J. Hodson, Diabetes 69, e10 (2020).
  * Speier and Rupnik [2003] S. Speier and M. Rupnik, Pflügers Arch. 446, 553 (2003).
  * Podobnik  _et al._ [2020a] B. Podobnik, D. Korošak, M. S. Klemen, A. Stožer, J. Dolenšek, M. S. Rupnik, P. C. Ivanov, P. Holme, and M. Jusup, $\beta$-cells operate collectively to help maintain glucose homeostasis: dataset (2020a), Archived at https://doi.org/10.17605/OSF.IO/NA5H3.
  * Podobnik  _et al._ [2020b] B. Podobnik, D. Korošak, M. S. Klemen, A. Stožer, J. Dolenšek, M. S. Rupnik, P. C. Ivanov, P. Holme, and M. Jusup, Biophys. J. 118, 2588 (2020b).
  * Slak Rupnik and Korošak [2019] M. Slak Rupnik and D. Korošak, Front. Physiol. 10, 1194 (2019).
  * Plerou  _et al._ [2002] V. Plerou, P. Gopikrishnan, B. Rosenow, L. A. N. Amaral, T. Guhr, and H. E. Stanley, Phys. Rev. E 65, 066126 (2002).
  * Johnson  _et al._ [2010] S. Johnson, J. J. Torres, J. Marro, and M. A. Munoz, Phys. Rev. Lett 104, 108702 (2010).
  * Newman [2002] M. E. Newman, Phys. Rev. Lett. 89, 208701 (2002).
  * Qu  _et al._ [2015] J. Qu, S.-J. Wang, M. Jusup, and Z. Wang, Sci. Rep. 5, 15450 (2015).
  * Albert  _et al._ [2000] R. Albert, H. Jeong, and A.-L. Barabási, Nature 406, 378 (2000).
  * Molloy and Reed [1995] M. Molloy and B. Reed, Random Struct. Algorithms 6, 161 (1995).
  * Nadakuditi and Newman [2013] R. R. Nadakuditi and M. E. Newman, Phys. Rev. E 87, 012803 (2013).
  * Majdandzic  _et al._ [2014] A. Majdandzic, B. Podobnik, S. V. Buldyrev, D. Y. Kenett, S. Havlin, and H. E. Stanley, Nat. Phys. 10, 34 (2014).
  * Podobnik  _et al._ [2017] B. Podobnik, M. Jusup, Z. Tiganj, W.-X. Wang, J. M. Buldú, and H. E. Stanley, Proc. Natl. Acad. Sci. USA 114, 11826 (2017).
  * Sun  _et al._ [2013] B. Sun, G. Duclos, and H. A. Stone, Phys. Rev. Lett. 110, 158103 (2013).
  * Potter  _et al._ [2016] G. D. Potter, T. A. Byrd, A. Mugler, and B. Sun, Proc. Natl. Acad. Sci. USA 113, 10334 (2016).
  * Fancher and Mugler [2017] S. Fancher and A. Mugler, Phys. Rev. Lett. 118, 078101 (2017).
  * Farkas  _et al._ [2001] I. J. Farkas, I. Derényi, A.-L. Barabási, and T. Vicsek, Phys. Rev. E 64, 026704 (2001).