Close
About
FAQ
Home
Collections
Login
USC Login
Register
0
Selected
Invert selection
Deselect all
Deselect all
Click here to refresh results
Click here to refresh results
USC
/
Digital Library
/
University of Southern California Dissertations and Theses
/
Understanding dynamics of cyber-physical systems: mathematical models, control algorithms and hardware incarnations
(USC Thesis Other)
Understanding dynamics of cyber-physical systems: mathematical models, control algorithms and hardware incarnations
PDF
Download
Share
Open document
Flip pages
Contact Us
Contact Us
Copy asset link
Request this asset
Transcript (if available)
Content
Understanding Dynamics of Cyber-Physical Systems: Mathematical Models, Control
Algorithms and Hardware Incarnations
by
Mahboobeh Ghorbani
A Dissertation Presented to the
FACULTY OF THE USC GRADUATE SCHOOL
UNIVERSITY OF SOUTHERN CALIFORNIA
In Partial Fulfillment of the
Requirements for the Degree
DOCTOR OF PHILOSOPHY
(ELECTRICAL ENGINEERING)
August 2019
Copyright 2019 Mahboobeh Ghorbani
Acknowledgements
First, I would like to express my gratitude to my thesis advisor, Paul Bogdan for his
guidance and advisement throughout the course of my Ph.D. work at the University of
Southern California.
Secondly, I consider myself fortunate to have received feedback on my research
from my thesis committee members: Professor Edmond Jonckheere, Professor Aiichiro
Nakano, Professor Sandeep Gupta. I am very thankful to Professor Aiichiro Nakano
for his perspective in Physics and suggestions about the implication of my gene study
work. I am very grateful for all the feedback and discussion related to the multi-variable
multi-fractal analysis and his feedbacks for gene network study. I am very thankful for
professor Sandeep Gupta for his suggestions on modeling model glucose dynamics. I
am also grateful for my qualifying exam committee members. Professor Bhaskar Krish-
namacharis challenging questions on the benefit of multi-variable multi-fractal analysis
helped me to look for more exciting application of the method on deciphering sleep
stages. I am also thankful for professor Massoud Pedram for his critical questions on
hardware realization of the artificial pancreas controller and his collaboration on cloud
computing workload analysis and control. I have been privileged to benefit from the
collaborations with Professor Yanzhi Wang and Professor Xue Lin at Northeastern Uni-
versity (when they were students at USC).
ii
My life as a graduate student at the University of Southern California would have
been much harder without the sincere support of many friends. I want to extend my
warm appreciation to the colleagues and friends who have made my Ph.D. life more
fruitful and exciting. Lastly, but not least, I want to present my heartfelt gratitude to all
the USC staff members (Diane, Liszl, Ted, Tim, Jeniffer, Tracy, Estela, Annie) who have
ensured that I have all the necessary support and advisement during my Ph.D. studies.
Above all, I want to share my deepest thoughts with my parents for their uncondi-
tional love and support. I am thankful to my mother for her emotional support, especially
and to my father for his thoughtful guides.
I am incredibly blessed to have the most profound love from my husband, my best
friend, Ehsan. I am thankful for his patience and encouragements during every step
of my study. Finally, I would like to express my gratitude to the agencies who have
contributed to the funding of my research, namely the National Science Foundation and
the University of Southern California.
iii
Contents
Acknowledgements ii
List of Figures vii
Abstract xii
1 Introduction 1
1.1 Medical Cyber-Physical Systems (CPS): definitions and challenges . . . 1
1.2 Theoretical foundations for CPS (fractals dynamics modeling) . . . . . 3
1.3 Thesis Overview: Research Objectives, Proposed Solutions and Future
Directions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
1.4 Thesis organization . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
2 Mathematical modeling and control algorithms for artificial pancreas 11
2.1 Mathematical modeling of blood glucose dynamics . . . . . . . . . . . 11
2.1.1 Statistical-Analysis-of-Blood-Glucose-Levels . . . . . . . . . . 11
2.1.2 Fractal Modeling and Optimal Control of Blood Glucose . . . . 14
2.2 Mathematical formulation of the control of blood glucose dynamics . . 17
2.2.1 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
2.2.2 Control related challenges and constraints . . . . . . . . . . . . 19
2.2.3 Control algorithms for BG level regulation . . . . . . . . . . . 21
2.3 Time-dependent fractal optimal control algorithm . . . . . . . . . . . . 27
2.4 Experimental results and discussion . . . . . . . . . . . . . . . . . . . 29
2.4.1 Model Identification and Analysis . . . . . . . . . . . . . . . . 29
2.4.2 Control Performance . . . . . . . . . . . . . . . . . . . . . . . 31
3 Comprehensive dynamic analysis of gene expression 34
3.1 Gene expression: definitions, modeling, and challenges . . . . . . . . . 34
3.2 Mathematical definitions of cross-dependence fractal dynamics . . . . . 36
3.2.1 The Hurst exponent and Multifractal Detrended Fluctuation Anal-
ysis (MFDFA) . . . . . . . . . . . . . . . . . . . . . . . . . . 36
3.2.2 Detrended Cross-Correlation Analysis (DCCA) . . . . . . . . . 37
iv
3.2.3 The Mandelbrot Binomial Cascade Model . . . . . . . . . . . . 38
3.2.4 Random Cascades on Wavelet dyadic trees . . . . . . . . . . . 39
3.2.5 Entropy and Entropy of a Network . . . . . . . . . . . . . . . . 40
3.3 Mathematical definitions of cross-dependence fractal dynamics . . . . . 41
3.3.1 Gene expression dynamics exhibits long-range dependency and
multifractal properties . . . . . . . . . . . . . . . . . . . . . . 41
3.3.2 The distribution of cross-correlation exponents of pairs of genes
and transcription factors of gene regulatory networks has a wide
range of variation . . . . . . . . . . . . . . . . . . . . . . . . . 44
3.3.3 Multifractal characteristics of interactions within the gene regu-
latory network can be modeled by random cascades on wavelet
dyadic trees . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44
3.4 Implications of cross-dependence fractal behavior in gene expression . . 46
4 Mutli-fractal multi variable detrended fluctuation analysis 57
4.1 Medical Complex (biological) systems: definitions and challenges . . . 57
4.2 Mathematical description of multi-fractal multi-variable detrended fluc-
tuation analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59
4.2.1 Description of the method . . . . . . . . . . . . . . . . . . . . 59
4.2.2 Computation complexity of mMFDFA method . . . . . . . . . 61
4.2.3 The correctness of the proposed Algorithm . . . . . . . . . . . 61
4.3 Verification of multi-fractal multi-variable detrended fluctuation analy-
sis on synthetic case studies . . . . . . . . . . . . . . . . . . . . . . . . 62
4.3.1 Application of the algorithm on independent time-series . . . . 63
4.3.2 Application of the algorithm on correlated time-series . . . . . 64
4.4 Verification of multi-fractal multi-variable detrended fluctuation analy-
sis on brain datasets . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68
4.4.1 Even a sleeping brain exhibits a complex coupled multi-fractal
dynamics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70
4.4.2 The multi-variable multi-fractal analysis reveals the distinction
of dynamics of deep sleep stages . . . . . . . . . . . . . . . . . 70
4.5 Verification of multi-fractal multi-variable detrended fluctuation analy-
sis on river network datasets . . . . . . . . . . . . . . . . . . . . . . . 71
4.5.1 Water flow dynamics across large river networks exhibits multi-
fractal characteristics. . . . . . . . . . . . . . . . . . . . . . . 71
4.5.2 Multi-fractality of water flow dynamics through a river basin is
season-dependent. . . . . . . . . . . . . . . . . . . . . . . . . 72
4.5.3 Multi-fractal spectrum of the river network dynamics encom-
passes the effect of climate change. . . . . . . . . . . . . . . . 74
5 Hardware prototyping of CPS components 76
5.1 Hardware design challenges in MCPS . . . . . . . . . . . . . . . . . . 76
v
5.2 Hardware prototyping of fractal optimal controllers . . . . . . . . . . . 77
5.3 Findings and suggestions for future work . . . . . . . . . . . . . . . . 80
6 Conclusions and future research directions 81
6.1 Major contributions of this Dissertation . . . . . . . . . . . . . . . . . 81
6.2 Limitations and future research directions . . . . . . . . . . . . . . . . 84
6.2.1 Complexity Aware Control . . . . . . . . . . . . . . . . . . . . 85
6.2.2 Controllability of Complex Networks . . . . . . . . . . . . . . 87
Reference List 88
vi
List of Figures
1.1 Fractional derivative as a non-local operator[1]. (a) First order deriva-
tive (b) fractional order derivative of order 0.8 (right). (c) log-log plot
of weights of W binomial coefficients (w
j
()) as a function of j. As can
be seen there is a power-law relation betweenw
j
andj. . . . . . . . . 5
1.2 Cyber Physical Systems: from sensing devices, to Modeling and Char-
acteristics extraction and Controlling algorithms . . . . . . . . . . . . . 7
2.1 (a). Inter-day (top panel), intra-day (middle panel) and inter-hour (bot-
tom panel) blood glucose dynamics for a patient suffering from type
1 diabetes mellitus (T1DM). Blood glucose time series were obtained
from DirecNet database [2]. Although inter- and intra-day blood glu-
cose variability is expected as a function of meal content and size, meal
absorption, level of physical and mental involvement in daily activities,
one can also observe a significant variation over a period of few hours
(bottom panel). (b). Inter-day (top panel), intra-day (middle panel) and
inter-hour (bottom panel) blood glucose variability for the person suffer-
ing from T1DM. One can clearly observe that the blood glucose incre-
ments display not only rich variability over a wide range of scales (from
days to hours), but also some degree of self-similar behavior. (c). Statis-
tical analysis of mean (top panel), variance (middle panel), and kurtosis
(bottom panel) shows that blood glucose dynamics cannot be modeled
adequately by linear state space models. . . . . . . . . . . . . . . . . . 12
2.2 (a). Empirical probability that the magnitude (absolute value) of the
negative blood glucose increments exceed a specific threshold deviates
from the Gaussian distribution and is better fitted by anstable distri-
bution. (b). Empirical probability that the positive blood glucose incre-
ments exceed a given threshold does not follow a Gaussian law and is
approximately fitted by anstable distribution. (c). The existence of
a fractal behavior and so inadequacy of modeling blood sugar dynam-
ics via linear state space models is proved by the multifractal spectrum
computed from the blood glucose time series of four individuals. . . . . 13
2.3 Systematic view of artificial pancreas [3]. . . . . . . . . . . . . . . . . 18
vii
2.4 Comparison between the blood glucose measurements and the time series
obtained from our proposed mathematical model for three cases: mag-
nification over a hypoglicemic period (a), approximately 900 blood glu-
cose samples (b), and magnification over a hyperglicemic period (c).
Our proposed mathematical model is able to capture well both the trend
and variability of blood glucose. . . . . . . . . . . . . . . . . . . . . . 29
2.5 a) Measured blood glucose level inmg=dL. b) Estimation error of best
ARMA model up to order 16. c) Estimation error of the proposed time
dependent fractal model. . . . . . . . . . . . . . . . . . . . . . . . . . 30
2.6 a) Root mean square error (RMSE) measured between the actual blood
glucose observations and an auto-regressive moving average (ARMA)
model. b) Root mean square error (RMSE) measured between the actual
blood glucose observations and the proposed time dependent fractal
model. c) Parameters , a, and b for the fractal model estimated over
non-overlapping time windows. . . . . . . . . . . . . . . . . . . . . . . 30
2.7 a) In-silico blood glucose level by applying integer and fractal optimal
control signals. b) Histogram of blood glucose levels of the integer order
controller due to Gaussian perturbations. c) Histogram of blood glucose
levels for the fractal controller in the presence of Gaussian perturbations. 33
3.1 (a). Gene expression process is initiated with the triggering of a gene in
the DNA strand and continues with the generation of RNA, mRNA and,
finally, the protein product. (b). Part of the gene regulatory network of
E. coli reported by [4], which consists of genes and transcription factors.
The diamonds represent transcription factors and the circles represent
genes. (c). Time series of transcription factors of E. coli obtained from
[5] and explained in the methods section. (d). Time series of genes of
E. coli obtained from [5] and explained in the methods section. . . . . . 52
3.2 Extracting a weighted graph from cross-correlation exponents of the
linked nodes in a gene regulatory network . . . . . . . . . . . . . . . . 53
3.3 (a). Scaling of fluctuation function of a gene time series in S. cerevisiae,
(b). Scaling of fluctuation function of ynel gene time series in E. coli,
(c). Histogram of cross-correlation exponent of genes in S. cerevisiae,
(d). Histogram of cross-correlation exponent of transcription factors in
S. cerevisiae, (e). Histogram of cross-correlation exponent of genes in
E. coli, (f). Histogram of cross-correlation exponent of transcription
factors in E. coli. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53
viii
3.4 (a). Scaling of detrended covariance function in a gene-transcription
factor link in E. coli (ihfB to ompR) and S. cerevisiae (YLR256W to
YKL020C) regulatory network, (b). Generalized cross-correlation expo-
nent of a gene-transcription factor link in E. coli and S. cerevisiae, (c).
Histogram of cross-correlation exponent of genes in S. cerevisiae, (d).
Histogram of cross-correlation exponents of gene regulatory network
links. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54
3.5 (a). Multifractal spectrum of several randomly picked links of gene
regulatory network of E. coli, (b). The multifractal spectrum of several
randomly picked links of gene regulatory network of S. cerevisiae, (c).
The multifractal spectrum of a link between the gene regulatory network
of S. cerevisiae and the best spectrum from the Mandelbrot binomial
model. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55
3.6 (a). The multifractal spectrum of a cross-dependency (link) between
the gene regulatory network of E. coli (soxR and soxS genes) that can
best be fitted by the log-normal cascade multifractal model, (b). The
multifractal spectrum of a cross-dependency (link) between the gene
regulatory network of S. cerevisiae (YKL032C and YKL043W genes)
that can best be fitted by the log-normal cascade multifractal model. . . 55
3.7 (a). The second derivative of the mass exponent for cross-dependencies
in E. coli. (b). The second derivative of the mass exponent for cross-
dependencies in S. cerevisiae. . . . . . . . . . . . . . . . . . . . . . . . 56
4.1 Simulation Results on independent time series:(a). Detrended fluctu-
ation function of the independent bi-variable signal along with detrended
fluctuation analysis of individual time series (b). The multifractal spec-
trum of independent time series (we have shown for bi-variable, 10, 40,
60 and 100 independent variables) . . . . . . . . . . . . . . . . . . . . 64
4.2 Simulation Results on correlated time series: (a). Detrended fluctu-
ation function of the correlated bi-variable signal along with detrended
fluctuation analysis of individual time series (b). The multifractal spec-
trum of correlated time series (we have shown for bi-variable, 10, 40, 60
and 100 correlated variables) . . . . . . . . . . . . . . . . . . . . . . . 66
4.3 The multifractal spectrum of correlated time series by different meth-
ods: a). The multifractal spectrum of correlated time series (we have
shown for a bi-variable case, 10, 40, 60 and 100 independent variables)
obtained by our proposed MMFDFA, (b). Multifractal spectrum of cor-
related time series obtained by the previous method proposed in [6, 7] . 67
ix
4.4 (a.) The multifractal spectrum of EEG time series for a healthy person
(n1) in S3 deep sleep state (b). The multifractal spectrum of EEG time
series for a healthy person (n1) in the REM sleep state (c). The mul-
tifractal spectrum of EEG time series for a healthy person (n2) in S3
deep sleep state (d). The multifractal spectrum of EEG time series for a
healthy person (n2) in a REM sleep state . . . . . . . . . . . . . . . . . 69
4.5 a) Map of the California water resources including Sacramento River
basins. b) Map of Minnesota water resources including upper Missis-
sippi River basins. c, d) Time series of the water flow for three differ-
ent measurement locations from the Mississippi (c) and Sacramento (d)
river basins in Oct 2007. e, f) Generalized Holder exponent of individ-
ual river flows and combination of individual time series in part (c) and
(d) respectively. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73
4.6 a) The generalized Hurst exponent as a function of the q-th order moment
for multi-channel river stream data obtained from four different sites
from the Mississippi river basin during the month of September in dif-
ferent years, i.e., 2012, 2013 and 2014 respectively. b) The general-
ized Hurst exponent (encompassing the degree of multi-fractality) as a
function of the q-th order moment for multi-channel river stream data
obtained from four different sites within the Sacramento river in Cali-
fornia during three different months in 2014. c) The generalized Hurst
exponent for the multi-channel river stream data obtained from four dif-
ferent sites of Sacramento River in California during the June-December
time period in several consecutive years, i.e., 2012, 2013, and 2014. . . 74
4.7 (a) And (b) Two sets of three time series representing the river flow
levels measured in the Mississippi river basin near Royalton, MN (Ch1),
at St. Cloud, MN (Ch2) and Mississippi River at HWY 610 in Brooklyn
Park, MN (Ch3) (c) And (d). The Generalized Hurst exponent as a
function of the q-th order moment for the single- and joint three-variable
time series data obtained from the Mississippi River locations shown in
part (a) and (b) respectively. . . . . . . . . . . . . . . . . . . . . . . . 75
5.1 (a). Primal-dual formulation of interior point method. (b). Pseudo-
code of LU factorization. (c). Interior point algorithm. (d). Parallel
architecture for updating every matrix element. . . . . . . . . . . . . . 78
5.2 (a). Clock frequency and (b). number of resources on FPGA for the
synthesized controller and different problem sizes. . . . . . . . . . . . . 79
x
6.1 a. Singularity spectra of the of heart rate signals for healthy and diseased
hearts ([8]). b. Complexity index (sum of Sample Entropy values for
scales 1 to 6, inclusively) and DFA-exponent (time scale greater than
one hour) for the non-diabetic and diabetic groups. Symbols with error
bars represent mean and SD values, respectively([9]). . . . . . . . . . . 84
6.2 Our proposed model predictive fractal control formulation. We have
incorporate the fractal dynamics into constraints of the problem in con-
strast to classical non-fractal optimal control problem. . . . . . . . . . . 86
xi
Abstract
Elucidating the dynamic nature of complex systems has significant consequences not
only on deciphering the topological interaction among various system components but
also for predicting their function and behavior. To date, most of the mathematical mod-
els for complex systems have ignored the spatio-temporal (multi-scale) correlations and
modeled their dynamics through memoryless (integer order differential or difference)
equations. Directly speaking, the current mathematical models of complex systems
neglect either the observed fractality in the dynamics of individual components of the
complex system or the higher-order cross-correlations between various components of
the complex system. To remedy this disconnect between the current mathematical mod-
eling and the intrinsic characteristics of complex systems, in this thesis, we address the
modeling of the dynamics and then, control of complex systems through fractal dimen-
sion analysis and fractional order control. To investigate the impact and importance of
multi-fractal characteristics in CPS design, we consider the design of the control algo-
rithm for artificial pancreas as a case study. We prove that the blood glucose (BG)
dynamics are highly complex and multi-fractal. We proposed a fractal-based model
for BG in contrast to previous memory-less models, and we show the advantage of our
model in terms of the goodness of fit. We formulate the fractional order model predictive
controller in the artificial pancreas in contrast to conventional integer order controllers.
xii
We have reported the advantage of the use of this model in bringing the BG level to
a reference point even in the presence of uncertainty. We have also shown the feasibility
of implementation of the fractal model predictive controller in hardware, and we have
reported the hardware prototyping results. We also extended our studies on gene expres-
sion networks as an example of a (biological) complex network. We investigate the scal-
ing properties of gene expression time series in E. coli and S. cerevisiae. We investigate
the individual gene expression dynamics as well as the cross-dependency between them
in the context of gene regulatory network. Our results demonstrate that the gene expres-
sion time series display fractal and long-range dependence characteristics. Besides, the
dynamics between genes and linked transcription factors in gene regulatory networks
are also fractal and long-range cross-correlated. The cross-correlation exponents in gene
regulatory networks are not unique. The distribution of the cross-correlation exponents
of gene regulatory networks for several types of cells can be interpreted as a measure of
the complexity of their functional behavior. We extend the concept of cross-correlation
fractal analysis and studied the dependencies through the concept of network entropy by
statistical distribution analysis on cross-dependencies. We also analyze the dynamical
process in complex networks. We investigate the existence of fractional order dynam-
ics in complex networks. We introduce a multidimensional (-variable) multi-fractal
detrended fluctuation analysis (mMFDFa) method that allows quantifying the correla-
tion complexity and memory exhibited between groups of interacting components of a
complex system. Our proposed method analyzes the higher-order statistics and quan-
tifies the power law scaling exhibited by dynamical fluctuations of coupled processes
(time series encompassing the evolution of complex system components). To illustrate
the depth and benefits of this mathematical method, we investigate the multi-fractal
spatio-temporal characteristics. First, we show how using this method is beneficial for
synthetic correlated and uncorrelated fractal signals. Then, we used this analysis to
xiii
distinguish the activity of human brain consisting of groups of neurons during various
sleep periods. As future direction toward this thesis, we propose to address a few new
problems. First, we want to formulate the complexity-aware control algorithm. Moti-
vated by the observed complex multi-fractal dynamics in physiological signals, there
can be a new class of control algorithms that can maintain the complexity of the target
control variable(s). Second, addressing the problem of controllability of complex net-
works considering fractal dynamics will also have a significant contribution to the field.
Lastly, investigating the scalability of the proposed fractal controllers and their hardware
implementation is also crucial to make the theory into practice.
xiv
Chapter 1
Introduction
1.1 Medical Cyber-Physical Systems (CPS): definitions
and challenges
Advances in embedded system research enables tight integration of several sensors mon-
itoring real world processes (e.g. heart rate, cloud movement, temperature fluctuation,
wind speed, earth magnetism fluctuations) and actuators able to control the environment.
These newly envisioned systems, expected to successfully integrate computation, com-
munication and control with physical systems are referred to as cyber physical systems
(CPS) [1] [2] [3]. The CPS denotes a network of embedded computational devices and
an associated set of wired or wireless networks than can monitor and control several
physical processes that occur in the environment. In the CPS area the goal is not only
to establish a reliable communication infrastructure between computational elements,
but also to include time and feedback-based control. Hence, a more direct interac-
tion between the system and physical world becomes possible. For example, vehicular
networks describing the cars movement in a city or the swarms of bacteria used for
diagnostic or drug delivery purposes are CPS examples that are distinct from classical
networked embedded systems.
Building a reliable CPS requires a new science of characterizing and control-
ling dynamic processes across heterogeneous networks of sensors and computational
devices. This science is needed to be able to model real world physical processes in
1
order to be able to reliably control them. Accurate characterization of workload proper-
ties (e.g. self-similarity and fractality) is a big challenge for CPS optimization.
Among all CPS domain, health care CPS need special attention. Cost of healthcare
in the U.S. is the highest in the world (per capita, GDP) [4]. Moreover, according to [2],
health costs across all OECD countries have increased and an important cause of this
is the use of more expensive and intensive treatments resulting from misuse, overuse
and/or underuse of care. Poor prevention, unhealthy life styles and lack of effective care
has led to higher rates of chronic diseases and related complications. For instance, in the
United States approximately 75 cents of every dollar is spent on patients with chronic
diseases. In 2005, this amounted to $1.5 trillion of the $2 trillion spent on health care.
Consequently, there is an urgent need for a better use of information technology for
monitoring, assessing and predicting the health state of people that would maximize
the success of prevention and treatment improving quality of life while minimizing the
health care costs related to hospitalization and side effects of medical care.
The use of information technology is already contributing in significant ways to
enhancing healthcare delivery and improving the quality of life. However, IT deploy-
ments have only scratched the surface of possibilities for the potential impact of com-
puter and information science and engineering on the quality and cost-effectiveness of
healthcare. Besides, with a vast amount of health-related data available through new
technologies and medical discoveries, there is massive potential for developing new
tools and approaches that will contribute to more cost-effective and sustainable health-
care. To address the increasing need for better biomedical systems, cyber physical
systems (CPS) represent the information technology quest of 21st century for a bet-
ter, cleaner, safer life, which integrates computation, communication and control with
physical resources. The CPS paradigm is meant to drive the design of many technolog-
ical systems ranging from smart transportation, to distributed energy and health care.
2
The health care CPS needs to maximize the patience’s quality-of-life, while minimizing
the intrinsic costs related to patient hospitalization.
1.2 Theoretical foundations for CPS (fractals dynamics
modeling)
The origin of fractional integral and derivative is almost as old as classical calculus itself.
In early 1695, in correspondence between L’Hopital and Leibniz in which they debated
the meaning of a 0.5 order differentiation [42], Leibniz wrote: thus it follows thatd
1
=2x
will be equal toL, ... from which one-day useful consequences will be drawn. It took
about 300 years for fractional calculus to be accepted as a practical instrument in physics
(e.g., visco-elasticity [34], dielectric polarization [39], heat transfer phenomena [6]) and
engineering (e.g., control [36], bio-engineering [33], traffic modeling [8], game theory
[11]).
Fractional (or fractal) calculus generalizes the notion of integer order derivative or
integral by dealing with functional differentiation and integration of arbitrary orders.
The first derivative of a function is approximated by the difference y=x where x =
h and y = f(x)f(xh). The second derivative is approximated by
2
y=
2
x,
where:
2
y = [f(x)f(xh)] (1.1)
= [f(x)f(xh)] [f(xh)f(x 2h)] (1.2)
=f(x) 2f(xh) +f(x 2h) (1.3)
Continuing this process, shows thatnth order derivative is approximated by :
3
d
n
f
dx
n
=
1
h
n
n
X
j=1
n
j
(1)
j
f(xjh) (1.4)
Where the binomial coefficients terms is defined as:
n
j
=
n!
j!(nj)!
(1.5)
The Grunwald definition of fractional derivative is the general form of equation 1.4,
where can take fractional values:
d
f
dx
=
1
h
X
j=1
j
(1)
j
f(xjh) (1.6)
Where the binomial coefficients are:
n
j
=
(n)
(j)(nj)
(1.7)
One interesting property of fractional calculus is that it is a non-local operator. It
means at each time x the fractional derivative depends on function values that are far
away from x as shwon in Figure 1.1.b. As can be seen, fractional derivative depends
on value of the function that have far distance in spite of the integer derivative (shown
in Figure 1.1.a which only depends on close neighborhood of the points. The contribu-
tion of value of the x
jh
point is measure by Grunwald weight w
j
= (1)
j
(n)
(j)(nj)
.
Figure 1.1.c shows the log-log plot of the w
j
with respect to j for different value
4
Figure 1.1: Fractional derivative as a non-local operator[1]. (a) First order derivative
(b) fractional order derivative of order 0.8 (right). (c) log-log plot of weights of W
binomial coefficients (w
j
()) as a function of j. As can be seen there is a power-law
relation betweenw
j
andj.
of order of fractional derivative namely . As depicted in the picture, the w
j
has a
power law decay function. This observation can also be verified by Sterling approxi-
mation [10] formula for factorial and gamma function. From Sterling approximation,
w
j
j
(1+)
().Hence the Grunwald definition is essentially a discrete convolution
with a power law:
Fractals
In what follows, we will briefly review the mathematical underpinnings behind
multi- fractal theory and provide the mathematical link between fractal behavior and
fractional order differentiation. First, let us recall the definition of multifractality for a
random function:
Definition: A real value random functionx(t)2 R is said to possess multifractal
property if for any real interval 0 <= h
min
<= h <= h
max
, there exists a continuum
set of fractal dimensions D(h) such that the function x(t) has a Holder-like behavior,
i.e.,:
5
jx(t +)x(t)j'jj
h
(x ! 0) (1.8)
Simply speaking the above definition has the following implications: (i) From a
mathematical perspective, the random functionx(t) satisfies the Holder condition and
is characterized by a set of Holder exponents which characterize the roughness of the
curve. Also, the Holder exponents are correlated with the scaling parameters of the
power laws observed at various scales of the random function. (ii) From a phenomeno-
logical perspective, the current state of the random function x(t) is affected by past
events weighted by power law functions (whose scaling coefficients depend on the
Holder exponents). Alternatively, one can state that for each exhibited Holder exponent
there exist a power law function weighing the past events and influencing the current
and future states of the random function. Mathematically, this can be easily expressed
through fractional calculus concepts such as fractional order derivative.
1.3 Thesis Overview: Research Objectives, Proposed
Solutions and Future Directions
Cyber-physical systems (CPS) are engineered systems that depend upon the synergy of
computational and physical components (see Figure 1.2). To build such cyber physical
systems, it is essential to have an accurate understanding of underlying physical systems
and then build engineered systems by engineering principals as depicted in Figure 1.2.
Toward this end, the proposed research addresses the following major objectives:
Modeling and analysis of physiological time series: Studying the statistical prop-
erties and dynamical features of physiological time series enable us to better
6
Figure 1.2: Cyber Physical Systems: from sensing devices, to Modeling and Charac-
teristics extraction and Controlling algorithms
understand the dynamics of the system and hence better prediction and control
ability. Toward this end, we have studied the dynamics of blood glucose time
series as a case study for the design of the artificial pancreas as an example of
health care cyber-physical systems. Also, we have studied the dynamics of gene
expression time series in E. Coli and S. Cerevisiae.
Modeling the dynamics of interactions inside physiological systems: We have
studied the couplings between gene expression time series in the gene regulatory
network as a case study. We have studies both pair-wise and network-wise inter-
actions. The result of this study was a statistical analysis of these interactions and
also new insight for multi-variable fractal analysis and network analysis. We pro-
pose a new algorithm for multi-variable multi-fractal time series analysis. Also,
7
we propose to extend the notion of network entropy to take the fractality of inter-
actions into account.
Control algorithm for cyber-physical systems: Based on the observed dynam-
ics for an underlying physical network, development of a new control algorithm
becomes necessary. We have formulated the fractal model predictive controller
in the design of artificial pancreas as a case study. We have shown how incorpo-
rating fractal dynamics in the design of control algorithm leads to better control
performance.
Hardware realization of the proposed control algorithm: We have investigated the
feasibility of hardware realization of the proposed control algorithm and proto-
typed a fractal model, predictive controller.
In summary, we have contributed mostly to the modeling and design of control algo-
rithm for CPS (the two upper layers in pyramid shown in Figure 1.2). As future direc-
tion, we propose formulation of new class of control algorithm that take into account the
complexity of physiological signals into objective function of control algorithms. Also,
we propose studying the complexity of dynamics of complex networks more deeply.
More explanations are present in the last Chapter of this thesis.
1.4 Thesis organization
This dissertation focuses on modeling, control, and realization of cyber-physical sys-
tems and mainly focuses on health care. Unlike linear control theory, we account for the
processes characteristics like multi-fractal behavior and present novel modeling, con-
trol, and implementation algorithms. In the following, we provide a brief overview of
our contributions.
8
Chapter 2 discusses first the importance of blood glucose control in the artificial
pancreas to improve the quality of life for diabetic patients. We show statistical and
multi-fractal characteristics of blood glucose, and based on that, we construct a finite
horizon model predictive multi-fractal controller. We show the impact of this controller
over the conventional controller. In summary, our contribution lies in modeling and
control of blood glucose, considering the multi-fractal nature of blood glucose time
series.
In Chapter 3, we study the dynamics of gene expression data in the context of gene
regulatory network for E. Coli and S. Cerevisiae. Our findings show not only individ-
ual gene expression time series were showing the fractal property, but also their cross-
dependency had this property. We extend the notion of network entropy and reported
the collective fractal dependencies of gene expression time series in Gene regulatory
networks.
In Chapter 4, we first show our developed new method for calculating the multi-
fractal spectrum of multi-variable time series. We show verification of our method
analytically and through several examples. Also, we compare our method to previous
works and show how they fail in distinguishing independent vs. correlated time series.
We also use our method on case studies, including studying human brain signals during
sleep and show how our approach is useful in deciphering deep sleep stages. We also
suggest using our method for studying the dynamics of complex networks in a broader
sense.
Chapter 5 highlights the importance of hardware realization of the control algorithm
proposed in Chapter 2 for the regulation of blood glucose. We discuss that this is a
case study, and we need to study the feasibility of hardware realization of hardware for
cyber-physical systems in general. Then, we show our hardware prototypings main idea
and implementation details on an FPGA device.
9
Chapter 6 elaborates the limitations on the methods developed by thesis and pro-
poses new research directions toward advancing them. We discuss how incorporating
the complexity of multi-fractal nature is only considered in the modeling part in our
approach and how integrating it as the objective of the control algorithm makes problem
statements more solid. We also discuss other problem to consider, which includes the
controllability of complex networks.
10
Chapter 2
Mathematical modeling and control
algorithms for artificial pancreas
2.1 Mathematical modeling of blood glucose dynamics
2.1.1 Statistical-Analysis-of-Blood-Glucose-Levels
Figure 2.1.a shows the blood glucose dynamics over three different time scales: a few
days (top panel), one day (middle panel) and a few hours (bottom panel). As one can
easily notice from Figure 2.1.a, the blood glucose level exhibits both a pronounced time
dependency and an irregular behavior characterized by fluctuations at all three time
scales. The time dependency can be explained by the fact that not all meals are of the
same type/size (meal variability) or are metabolized by the organism in the same way
throughout the day (metabolic variability). Moreover, there are many other factors influ-
encing blood glucose dynamics such as stress and exercise levels (lifestyle variability).
An alternative perspective to the nonlinear dynamics of blood glucose is shown in Figure
2.1.b, where we magnify over the magnitude of fluctuations (i.e., blood glucose incre-
ments are computed as followsg(t)g(t 1)) over the same three time scales. This
erratic behavior observed in blood glucose fluctuations corresponds to the mathematical
concept of statistical self-similarity
1
.
1
The statistical self-similarity property of a stochastic process implies that its probability distribution
function (PDF) and all its statistical moments are invariant to mathematical scaling; simply speaking, sta-
tistical self-similarity implies that some degree of repetitions can be observed in the stochastic trajectories
over multiple time scales.
11
Figure 2.1: (a). Inter-day (top panel), intra-day (middle panel) and inter-hour (bot-
tom panel) blood glucose dynamics for a patient suffering from type 1 diabetes mel-
litus (T1DM). Blood glucose time series were obtained from DirecNet database [2].
Although inter- and intra-day blood glucose variability is expected as a function of meal
content and size, meal absorption, level of physical and mental involvement in daily
activities, one can also observe a significant variation over a period of few hours (bot-
tom panel). (b). Inter-day (top panel), intra-day (middle panel) and inter-hour (bottom
panel) blood glucose variability for the person suffering from T1DM. One can clearly
observe that the blood glucose increments display not only rich variability over a wide
range of scales (from days to hours), but also some degree of self-similar behavior. (c).
Statistical analysis of mean (top panel), variance (middle panel), and kurtosis (bottom
panel) shows that blood glucose dynamics cannot be modeled adequately by linear state
space models.
There are two implications of statistical self-similarity in blood glucose dynamics:
First, the blood glucose can have a complex behavior characterized by non-Gaussian
distribution and so the linear time invariant system theory is not adequate for modeling
and controlling its dynamics. Indeed, by analyzing the mean, variance and kurtosis of
blood glucose over time windows of five hours, we observe not only that blood glucose
level vary greatly over time, but also that blood glucose variability is not well approx-
imated by Gaussian distribution
2
. In addition to this, one can clearly see from Figures
2.2.a and 2.2.b that both the negative (in absolute value) and positive increments are
better fitted by the -stable distribution rather than by the Gaussian law. Second, the
2
A Gaussian stochastic process is characterized solely by the first two moments (i.e., mean and vari-
ance) and all its higher order moments are zero.
12
Figure 2.2: (a). Empirical probability that the magnitude (absolute value) of the negative
blood glucose increments exceed a specific threshold deviates from the Gaussian distri-
bution and is better fitted by anstable distribution. (b). Empirical probability that
the positive blood glucose increments exceed a given threshold does not follow a Gaus-
sian law and is approximately fitted by anstable distribution. (c). The existence of
a fractal behavior and so inadequacy of modeling blood sugar dynamics via linear state
space models is proved by the multifractal spectrum computed from the blood glucose
time series of four individuals.
exhibited self-similar behavior of blood glucose may be indicative of a fractal behavior.
By applying the box counting method
3
in an iterative manner over windows of blood
glucose time series, Figure 2.2.c shows the distribution of fractal dimensions. The wide
range of exhibited fractal dimensions shows both that blood glucose exhibits a much
more complex dynamics that cannot be approximated by linear state space equations
and that the mathematical models need to account for time dependent fractal behavior
(i.e., for higher accuracy, the mathematical modeling of blood glucose should incorpo-
rate a time dependent parameter representing the fractal behavior and so the level of
long range memory present in the organism as a result of various inter-related metabolic
processes). This will be discussed in next section.
3
Box counting method refers to a mathematical algorithm for analyzing complex patterns in space and
time and quantifying whether the considered object or stochastic process possesses fractal characteristics.
The box counting procedure investigates how the number of boxes of a certain size required to cover an
object vary as a function of the box size and determines the spectrum of fractal dimensions [11, 12].
13
2.1.2 Fractal Modeling and Optimal Control of Blood Glucose
Patients with type-1 diabetes have little or no insulin production, leaving the body
unable to lower the blood glucose levels without exogenous insulin. The AP purpose
is to rely on an accurate modeling of the blood glucose dynamics, determine the right
(best) insulin injection rate and infuse the insulin via a specialized pump such that the
blood glucose level is kept within a normal range. Toward this aim, in this section, we
first build a novel mathematical model able to capture the blood glucose characteris-
tics observed in real world measurements (see Figs. 2.1 and 2.2) via fractional calculus
[13] concepts. Building on our proposed time dependent fractal model of blood glucose
dynamics and using concepts from optimal control theory, we formulate a constrained
fractal optimal control problem which finds the best insulin injection rate that brings the
magnitude of blood sugar within the normal range of a healthy individual.
To date, the state-of-the-art in AP control algorithms relied on the assumption that
the auto-correlation (i.e., two-time correlation) function of the blood glucose dynamics
presents an exponential decay corresponding to a memoryless stochastic process and
ignored the intrinsic physiological variability. In addition, most current control algo-
rithms relied on memoryless linear models to find the amount of insulin necessary to
be injected for canceling the glucose rise due to metabolic variability. Nevertheless, the
quality and performance of such control algorithms are highly dependent on the accu-
racy of modeling blood glucose dynamics. Higher model accuracy, especially for non-
linear dynamics, is not synonymous with increased number of parameters in the model
(e.g., increasing the order of ARMA may decrease the modeling errors, but does not
allow for better accuracy and predictability [14]). Instead, by accurately characterizing
the complexity of a physiological process via nonlinear (higher order) temporal corre-
lations allows us to identify compact and robust nonlinear models which provide a high
degree of accuracy with minimum number of parameters. Hence, in what follows, we
14
propose a compact phenomenological model of blood glucose dynamics starting from
statistical characterization of exhibited glucose fluctuations.
The statistical analysis reported in Figs. 2.1 and 2.2 show that blood glucose level
exhibits a complex time dependent fractal behavior (i.e., blood glucose presents some
degree of long-range correlations which is in contrast to the short-range memory prop-
erty of current models). In other words, blood glucose fluctuations occurring on multiple
time scales are correlated over considerable periods of time (e.g., tens of hours and even
days) and so posses long range memory properties
4
. Consequently, rather than ignoring
these characteristics, we develop a phenomenological mathematical modeling of blood
glucose dynamics capable of encapsulating the time dependent fractal behavior observed
in the spectrum of fractal dimensions shown in Fig. 2.2.c.
To make the discussion more concrete, we define byg(t) (wheret > 0) the blood
glucose level at time t. To capture the time dependent fractal behavior shown in Fig.
2.2.c, and implicitly the long range memory property, we introduce the fractional deriva-
tive (i.e.,
0
D
(t)
t
) of the blood glucose levelg(t) and express its dynamics as follows:
0
D
(t)
t
g(t) =a(t)g(t) +b(t)u(t) +c(t)d(t) (2.1)
where a(t) is the , u(t) denotes the amount of insulin injected at time t, b(t) is a
coefficient representing the impact of injected insulin on the blood glucose dynamics,
c(t) and d(t) are error coefficient and error term assumed to affect the blood glucose
dynamics at timet, and(t) is the order of the fractional derivative introduced to model
the time dependent memory of the blood glucose dynamics and is defined as follows:
4
For long range dependent processes, the coupling between values at different times decreases as a
power law of the time differences.
15
0
D
(t)
t
g(t) =
1
(n(t))
d
dt
n
t
R
0
(t)
n(t)1
g()d (2.2)
withn being the smallest integer upper bounding(t), is the gamma function [15][16].
From a phenomenological perspective, the time dependent fractional differentiation
0
D
(t)
t
g(t) implies that the current and future blood glucose levels depend on a set of
previous values and so it captures the observed long range memory and fractal behavior
observed in the power law glucose fluctuations and the spectrum of fractal dimensions
shown in Fig. 2.2. In addition, rather then iteratively searching for the best order of an
ARMA model and so increasing the number of time dependent parameters in the model,
here we only introduce (t) to capture the effect of past values on current and future
blood glucose dynamics. Of note, the introduction of a time dependent fractal (fractional
order) exponent(t) offers also a compact representation of the exhibited intra- (due to
meal absorption, subcutaneous insulin absorption, daily physical exercise, stress) and
inter- (circadian rhythm, weekly activities) day variability. Consequently, the proposed
model in (2.1) is able to capture the blood glucose variability more adequately then
integer order differential equations and addresses the modeling challenges emphasized
by the Diabetes Technology Meeting [17].
16
2.2 Mathematical formulation of the control of blood
glucose dynamics
2.2.1 Motivation
In healthy individuals, the alpha and beta cells of the pancreas regulate the blood glu-
cose concentration to around 80 mg/dl. For people suffering from Type-1 diabetes mel-
litus (T1DM), which is one of the fastest growing diseases globally, there is little or no
endogenous insulin production, leaving the body unable to lower blood glucose with-
out exogenous insulin. The impact of the intensive insulin therapy was not revealed up
until the publication of results of Diabetes Control and Complication Trial (DCCT) in
1993 [18]. The CDDT involved a comparison of conventional therapy (one or two daily
insulin injections and a daily monitoring of blood glucose or urine) and intensive insulin
therapy and concluded that intensive therapy resulted in lower mean blood glucose
values and significantly reduced complications (retinopathy, nephropathy and macro-
vascular disease). The risk of complication is directly related to glycated hemoglobin
known as HbA1c. OGrady et al. find that tighter blood glucose levels achievable with
a closed-loop artificial pancreas (AP) results in Medicare savings of 1.9 billion over 25
years with improved quality of life (QOL) [19]. A schematic view of a closed-loop
artificial pancreas is shown in Figure 2.3, which is mainly composed of three parts:
Continuous time blood glucose measurement (CGM): The knowledge of glucose
concentration in blood is a key aspect in the quantitative understanding of the glucose-
insulin system and in diagnosis and treatment of diabetes. By the ability of CGM devices
to provide glucose readings in real time, engineers can exploit signal processing and con-
trol theory to be used in designing efficient artificial pancreas. Besides the improvement
of hardware part of CGM devices, a vast amount of research has devoted to address
denoising, prediction and alert generation [20][21][22][23][24][25][26].
17
Figure 2.3: Systematic view of artificial pancreas [3].
Control algorithm and safety layer: The main component of AP is the control
algorithm that determines the right insulin injection rate based on CGM data to pre-
vent hyperglycemia and hypoglycemia. Design of a QOL-aware AP is a challenging
task since it requires building accurate mathematical model of glucose-insulin kinetics.
Algorithms often include a safety layer as a supervisory module that constraints insulin
delivery. This layer may monitor and limit insulin on board (the insulin delivered but
yet to exert its action) or maximum insulin rate or may suspend insulin delivery at low
glucose levels or when glucose is decreasing rapidly. In this paper, we overview these
challenges and control strategies employed so far.
18
Insulin injection device: The essential function of AP is the Insulin delivery.
Insulin pumps, if inserted in a proper closed-loop system allow automatic insulin deliv-
ery. There are several technologies that can perform this task: an intra-venous route,
subcutaneous insulin infusion (SCII) or intaperitoneal insulin delivery. Continuous sub-
cutaneous insulin infusion (CSII) uses a portable electromechanical pump to mimic non-
diabetic insulin delivery as it infuses at preselected rates normally a slow basal rate with
patient-activated boosts at mealtime.
2.2.2 Control related challenges and constraints
Blood glucose (BG) regulation requires control algorithm to determine the best insulin
injection over time. They have been tested in-silico and clinically over time and
improved over the years. In this section, we first address the main challenges and con-
straints in designing efficient control algorithm. Next, we explain different control algo-
rithms proposed so far and compare their performance.
Non-negligible delay in glucose measurement and between insulin injection and
absorption: After administration of a subcutaneous bolus of rapid acting insulin ana-
logues, the maximum BG lowering effect may occur after up to 90-120 min. This time
lag is often not accounted for design of control algorithm. Patients treated with insulin
pump are warned against stacking caused by the administration of a series of correc-
tion boluses. The same principle applies to closed-loop systems. In order to prevent
hypoglycemia, high glucose levels have to be brought within normal range slowly dur-
ing closed-loop delivery. Methods to assess the impact of injecting insulin (e.g. the
one proposed in [27]) are highly needed in order to protect against insulin overdosing.
Two alternative insulin delivery routes, intraperitoneal (IP) and technosphere insulin
(TI) showed faster pharmacokinetic characteristics that can improve the design of future
19
AP systems. Design of the AP using these fast acting alternative routes may enhance
BG regulation by reducing actuation delays, especially during mealtime.
Asymmetric risk for low and high BG levels: The ultimate objective of any AP is
to improve QOL and minimize complications resulting from poor blood glucose control.
Toward this end, one should note to the asymmetric risk associated to high BG levels.
Low BG levels are acutely risky as they can result in altered mental state, seizures and
coma. Meanwhile, high BG levels increase the risk of chronic complications such as
retinopathy, nephropathy and cardiovascular disease.
Irreversible action of insulin: Only positive amount of the injected insulin is pos-
sible and it cannot be collected back from the patients blood. An alternative to deal with
this problem is to use bihormonal treatment [28] consisting of injecting glucagon and
insulin. However, this also increases the problem space and complexity.
Meal detection/estimation: Meal dynamics can have a significant disturbance
effect on BG level. In a fully closed-loop mode, insulin is delivered on the basis of
glucose excursions only, without information about timing or meal size. In a less ambi-
tious configuration that uses meal announcement, the closed-loop system is informed
about meal size, and may generate advice on prandial insulin bolus.
Alternatively, control algorithms can automatically increase insulin delivery based
on the carbohydrate content of the meal. A hybrid approach is characterized by admin-
istration of a small pre-meal priming bolus or administration of a fixed bolus and deliv-
ering the remaining insulin through the closed-loop operation [29].
Time dependency of control requirements: An important challenge in develop-
ment of artificial pancreas is that overnight treatment requires slow acting insulin injec-
tion while post-prandial control requires rapid and aggressive insulin delivery to control
BG.
20
On the other hand, exercise of moderate intensity increases the risk of hypoglycemia
[30]. Exercise announcement or heart rate monitoring to suspend insulin during closed-
loop delivery may be another effective method to control glucose levels during exer-
cise. Preemptive carbohydrate intake or dual hormone treatment with glucagon might
be needed to fully eliminate the risk of exercise-related hypoglycemia as responses to
exercise are highly variable. To sum up, BG control is a time dependent process and this
should be taken into account in order to have a safe and efficient AP.
Variability of model parameters: Up to 4 times inter-subject variability in rapid-
acting insulin analogue pharmacokinetics has suggested with occasionally as much as
50 % intra-subject variability [28]. Within subject variability of insulin needs includes
both day-to-day and hour-to-hour variations in insulin sensitivity owning to circadian
and diurnal cycles, dawn phenomenon (an abnormal early morning increase in BG con-
centration), acute illness, stress and a delayed effect of alcohol intake. Basal insulin
needs are generally lower in young individuals compared with older ones. Also, since
overnight control requires regulation based on mild control actions while postprandial
regulation is characterized by prompt and energetic correction, timely control effect
should take place.
2.2.3 Control algorithms for BG level regulation
In this section we present two main groups of controller for BG level regulation namely
proportional-integral-derivative and model predictive controller.
Proportional-Integral-Derivative (PID) controller
PID controller is a generic control loop feedback mechanism widely used in industrial
control. The PID control algorithm for artificial pancreas adjusts the insulin delivery rate
by assessing glucose excursions from three viewpoints: the departure from the target
glucose level (the proportional component), the area under the curve between measured
21
and target glucose levels (the integral component) and the rate of change in measured
glucose levels (the derivative component). Some controllers include only a subset of the
components (e.g. a proportional-derivative [31]).
To better understand the intuition behind using PID controller in the control algo-
rithm for artificial pancreas, one should note that dose of insulin is directly related to the
proportional error (P) (current glucose minus target glucose) since a patient with higher
glucose level needs more insulin rather than one with lower glucose level. Moreover,
in two patients with the same glucose level but with different rate of glucose increase,
the one with higher increase rate should get a higher dose of insulin and this justifies
the derivative element (D). To understand the role of Integral element (I), it should be
noted that for two patients with the same current glucose level and no change in the very
recent minutes, the one with more hours spending in high BG level (thus, having more
integral error) needs more insulin due to the fact that this is a sign on insulin resistivity.
Steil et al. have shown the normal healthy pancreas displays proportional, derivative
and integral dynamics [32]. They argue that the abrupt step increase in glucose causes a
rapid rise in pancreatic insulin release, which is called first phase response and is related
primarily to the derivative component. Slower rise in insulin is called the second phase
response, which corresponds to proportional term and persists as long as glucose is ele-
vated. There is also an integral component employed in the second phase since insulin
secretion after 3 hours of elevated glucose at a fixed level is greater than insulin secretion
after only 1 hour at the same glucose level.
Equation 1 shows the the components of control signal (u(t)) that is the amount
of insulin injection rate as a function ofe(t) which is the difference between BG level
and the reference value. The K
p
, K
i
and K
d
parameters can be assigned by learning
algorithms that have been discussed in control related textbooks. Optimizing using PID
22
controller needs tuning of the controller by some methods, like the ones proposed in
[33] and [34].
u(t) =K
p
e(t) +K
i
t
R
0
e()d +K
d
de(t)
dt
(2.3)
PID approach has inherent limitations due to time lags in glucose sensing and insulin
action. Several studies have investigated this approach and achieved some improvement
over conventional PID approach. For example, Weinzimer et al. in [35] have tested
PID algorithm for insulin injection in 17 adolescences. They have tested both fully
closed-loop and hybrid closed-loop (with pre-meal priming bolus) and show the addi-
tion of small manual priming bolus doses of insulin given 15 min before meals improves
glycemic excursions. A different study by Renard et al. in [36] proved the feasibil-
ity of intraperitoneal insulin delivery for artificial beta cell and supported the need for
further study since subcutaneous insulin delivery from a portable pump encountered
delays and variability in insulin absorption. They evaluated their proposed method in a
clinical study on eight T1DM patients while the time spent in 4.4-6.6 mmol/l was the
primary end point. Another study in [28] uses both insulin and glucagon to prevent
hypoglycemia encountered in PID algorithm with only insulin as the treatment.
Model Predictive Controller
Model Predictive Control (MPC) is a general optimization framework that can involve
many different types of models and objective functions. The MPC approach is at the
front of current research into closed-loop systems. It acceptably accommodates delays
associated with insulin absorption and can also easily account for meal intake and pran-
dial insulin boluses by the patient. The other advantage of model predictive control
paradigm is the fact that it can account for variability since the model parameters can be
personalized. The main advantage of MPC is the fact that it allows the current timeslot
23
to be optimized, while keeping future timeslots into account. This is achieved by opti-
mizing a finite time-horizon, but only implementing the current timeslot. MPC has the
ability to anticipate future events and can take control actions accordingly.
The vital ingredient of MPC is a model that links insulin delivery and meal ingestion
to glucose excursions. This model can be physiological and account for fundamental
processes regulating glucose levels or a black box model that disregards insights but
learns the insulin glucose relationship via formal pattern recognition technique. They
both can benefit from a wide range of mathematical models of glucoregulatory system.
It is therefore clear that proper models of glucose and insulin kinetics as well as models
that can be used to predict near-future metabolic behavior are mandatory. Minimal
models (describing the key components of system functionality) and maximal models
(nonlinear, high order models) are reviewed by Cobelli in [26].
A general MPC problem formulation, which includes optimization objective (Equa-
tion 2.4), glucose-insulin dynamical model (Equation 2.5) and initial value, glucose state
and insulin control constraints (Equation 2.6) can be written as follows.
min
u(t)
t
f
R
0
F (g(t);u(t))dt (2.4)
dg(t)
dt
=a
G
g(t) +b
G
u(t) (2.5)
g(t = 0) =g
0
;u
min
u(t)u
max
;g
min
g(t)g
max
(2.6)
whereg(t) denotes the BG level andu(t) denotes the amount of insulin injected at
timet which should be determined by solving the optimization problem;a
G
andb
G
are
coefficients representing the impact of injected insulin on the BG dynamics. Also, t
f
24
represents the finite horizon of the control problem which is usually 2h to 4h predic-
tion window that corresponds to the bulk duration of action of a rapid acting insulin
analogue such as aspart, lispro and glulisine. g
ref
(t) is the time dependent glucose ref-
erence value that can be chosen depending on the current state to avoid hypoglycemia or
hyperglycemia. Initial condition is addressed by includingg
0
which is the initial glucose
level. Finally,u
min
andu
max
are the minimum and maximum allowed insulin amounts
to be injected andg
min
andg
max
are the lower and upper bounds on the glucose level.
F (g(t) is a generic form for all possible cost functions. But, it is usually desired to
minimize a summation form including both the distance to the reference glucose value
and insulin injection effort. MPC has shown to be suitable for multivariate nonlinear
systems such as the human body and it significantly gives better performance than PID
control with patient-specific tuning. Several variations of MPC have been proposed in
the literature. We briefly categorize them as follows:
Linear model predictive control (LMPC): There are several research studies that use
linear model predictive controller. The work presented in [37] was the first in silico
trial for linear model predictive approach which also showed better performance of
MPC rather than PID controller in terms of limiting the oscillation of glucose levels.
Research efforts presented in [38] and [29] were the first clinical investigations of linear
model predictive algorithms in artificial pancreas that reported the superiority of using
this approach over PID controller. in [37] Magni et al. in present an unconstrained MPC
where the model is a linearization of a nonlinear parameters.
Researchers extend the MPC by defining new types of objective function and addi-
tional features to the problem formulation. Heusden et al. proposed using a priori patient
characterization and fitting a linear control relevant model around the control point [39].
They also defined a new cost function named zone model predictive in contrast to pre-
vious studies in which only the distance to a target reference point is considered as the
25
cost function. They consider a range for the BG level as the objective of the optimization
and define the cost function as the minimum distant to the preferred zone. They have
verified the robustness of the algorithm in silico and showed that the hypoglycemia is
completely avoided even after meal disturbance. Lee et al. in [40] use new meal size
estimation algorithm to the integrated AP and show how its performance is better than
MPC-only case.
Non-linear model predictive (NMPC): Hovorka et al. in [41] present a nonlinear
model and Bayesian techniques to estimate parameters in simulation studies. Clini-
cal studies were performed under fasting conditions based on measurements that were
delayed by 30 min to mimic the time lag associated with a sensor. The authors per-
formed overnight studies using an algorithm and transferring results to a pump at 15
min intervals. The major result was a reduction in nocturnal hypoglycemia compared
to standard pump treatment. Zarkogianni et al. also use a nonlinear model-predictive
control for prediction of BG and control algorithm [42]. They have shown the useful-
ness of using this nonlinear MPC in silico for different meal profiles, fasting conditions,
inter-patient variability and intraday variation.
Fractal model predictive control (FMPC): In spite of significant amount of work
in PID and MPC, the complexity of BG dynamics has not been fully addressed. For
instance BG is time dependent process that is influenced by various factors (meal size,
exercise, psychological state, etc.). This has prompted a comprehensive multifractal
investigation of BG dynamics [3] from publicly available data set [43]. The authors
have shown how using fractional order controller leads more robust control over con-
ventional integer order model predictive controller. They formulate the BG dynamics
as a time dependent fractional order control problem and report the feasibility of imple-
mentation of fractional controller in hardware and report their results in terms of area
26
and speed in field programmable gate array (FPGA). We compare the impact of apply-
ing fractional order controller to the conventional first order derivative controller. Fig
2 shows the outcome of applying both types of controllers to bring to some reference
value which is 100mg=dL in this case. Unlike the expectation of integer order controller
the final glucose value at the end of control horizon is much lower than the one expected.
2.3 Time-dependent fractal optimal control algorithm
The key part of the AP is to robustly and safely regulate the blood glucose by insulin
injection such that either hypoglcemia or hyperglicemia are avoided. Consequently,
we formulate the blood glucose regulation as a time dependent fractal optimal control
problem:
min
u(t)
t
f
R
0
(g(t)g
ref
(t))
2
+w (u(t)u(t 1))
2
dt (2.7)
subject to glucose-insulin dynamics model in (2.8), initial value (2.9), glucose state
(2.10) and insulin control (2.11) constraints:
0
D
(t)
t
g(t) =a(t)g(t) +b(t)u(t) +c(t)d(t) (2.8)
g(t = 0) =g
0
(2.9)
u
min
u(t)u
max
(2.10)
g
min
g(t)g
max
(2.11)
wheret
f
represents the finite horizon of the control problem,g
ref
(t) is the time depen-
dent glucose reference value,w is an insulin sensitivity coefficient that both weighs the
27
importance of the two mathematical terms and accounts for the different units of mea-
sure of the two terms,u(t) denotes the insulin injected at timet,g
0
is the initial glucose
level, u
min
and u
max
are the minimum and maximum allowed insulin amounts to be
injected,g
min
andg
max
are the safe bounds imposed on the glucose level.
The aim of the controller is to find the near future optimal insulin amounts which
minimizes both the deviations of the blood glucose level from the predefined reference
values and the amount of insulin injected at each timet. At each control interval only the
first insulin valueu is injected through the insulin pump and then the optimization prob-
lem is solved again to determine the next amount for the newly observed blood glucose
dynamics. Of note, we choose to minimize the squared difference between consecutive
insulin injection amounts (u(t)u(t 1))
2
rather than (u(t))
2
, because. Nevertheless,
our controller can easily be modified to minimize the deviations between the injected
insulin and a predefined clinical insulin injection profile (i.e., (u(t)u
ref
)
2
).
L =
t
f
R
0
f(g(t)g
ref
(t))
2
+w (u(t)u(t 1))
2
+
1
(t)
h
0
D
(t)
t
g(t)ag(t)bu(t)cd
i
+
2
(uu
min
) + (2.12)
3
(u
max
u) +
4
(gg
min
) +
5
(g
max
g)gdt
and imposing the following conditions:
@L
@g
+
t
D
(t)
t
f
@L
@
0
D
(t)
t
g
= 0;
@L
@u
= 0;
@L
@
k
= 0; k = 1 5 (2.13)
Towards this end, the control algorithm consists of discretizing the optimality con-
ditions summarized in (2.13) over the considered finite horizon [0;t
f
] and solve a linear
program that determines the right amount of insulin to be injected such that the devia-
tions from the reference profile are minimized.
28
Figure 2.4: Comparison between the blood glucose measurements and the time series
obtained from our proposed mathematical model for three cases: magnification over a
hypoglicemic period (a), approximately 900 blood glucose samples (b), and magnifi-
cation over a hyperglicemic period (c). Our proposed mathematical model is able to
capture well both the trend and variability of blood glucose.
2.4 Experimental results and discussion
2.4.1 Model Identification and Analysis
Building on the non-stationary and fractal nature of blood glucose dynamics, we devel-
oped an efficient wavelet-based parameters estimation method for the mathematical
model in eq. 2.1. Temporal resolution of the wavelet transform allows us to capture
short-lived time phenomena’s like singularity point. Also, since wavelets are localized
in scale they identify the self-similar (scaling) behavior in stochastic processes. Starting
from the theoretical premises of wavelets, the model identification consists of two steps:
29
Figure 2.5: a) Measured blood glucose level in mg=dL. b) Estimation error of best
ARMA model up to order 16. c) Estimation error of the proposed time dependent fractal
model.
Figure 2.6: a) Root mean square error (RMSE) measured between the actual blood
glucose observations and an auto-regressive moving average (ARMA) model. b) Root
mean square error (RMSE) measured between the actual blood glucose observations
and the proposed time dependent fractal model. c) Parameters,a, andb for the fractal
model estimated over non-overlapping time windows.
first, a wavelet scaling analysis inspired from [44] is used to estimated the time depen-
dent order of the fractional derivative in eq. (2.1); second, a low dimensional
5
linear
regression problem is solved efficiently to find the model parametersa,b andc.
5
Of note, by using the wavelet representation of the observed blood glucose measurements, we use
only the wavelet coefficients at each scale rather than each and every measurement. This reduces the size
of the linear regression by two orders of magnitude.
30
To validate the accuracy of our proposed model in (2.1), we used blood glucose
time series of 10 patients from publicly available databases [2]. We compared in Fig.
2.4 the accuracy between previous auto-regressive (AR) model (up to 16 order), an
improved auto-regressive moving average (ARMA) (up to 16 order), and our fractal
model (2.1) in terms of root mean square error (RMSE) and maximum deviation (MD)
measured between the observed data and extracted model. For illustration purposes,
Fig. 2.5 shows the comparison between time series obtained from our model and the
measured blood glucose levels in one of the patients. Also, we have plotted the estimated
parameters and the difference of the estimated fractional model from the real data and
the difference of the best ARMA model estimation in Fig 2.6.c. As can been seen from
the Figs. 2.6.a and 2.6.b, the proposed non-stationary fractal model is superior to the
conventional ARMA model. More important, the non-stationary model has consistently
smaller RMSE than ARMA and so is able to better capture the sudden changes in the
blood glucose. This represents a significant benefit for health-care purposes where the
anomaly detection is important for controller devices.
2.4.2 Control Performance
To illustrate the performance of the time dependent fractional order controller discussed
in Section 3.2, we consider the blood glucose abnormal dynamics for an individual suf-
fering from type-1 diabetes (real measurements from [2]). We apply the controller when
the patient’s blood glucose exhibits a highly variable as the one in Fig. 2.5.c and mea-
sures 150mg=dL which is above the maximum normal range of 140mg=dL. As the
normal range of the blood glucose for a healthy individual is between 80mg=dL and
140mg=dL, we set the time dependent reference value g
ref
(t) to be 100mg=dL. The
reason for choosing this reference value is two fold: i) First, we aim at stabilizing the
blood glucose level that is just slightly above the normal range of 80mg=dL such that
31
we avoid that upon control action due to uncertainty in insulin absorption the individual
will not experience a hypoglicemia. ii) Second, we also seek to lower the blood glu-
cose because of detrimental long term effects of hyperglicemia that is experienced for
prolonged times. For instance, increased level of blood glucose due to hyperglicemia
amplifies the heart tissue damages, provokes coronary endothelial dysfunction and raises
the likelihood of myocardial ischemia [45].
To compare the impact of the fractional order controller to the conventional first
order derivative controller, we have applied the fractional order and conventional opti-
mal control algorithms on the observed time series (see Fig. 2.7). First, we identified an
ARMA based model of the blood glucose dynamics and solved the conventional (integer
order) optimal control problem to find the non-fractal control signalu
nonfractal
(t). We
applied the newly computed control action u
nonfractal
(t) to the most accurate blood
glucose dynamics model and determined the stabilization level of the glucose as a result
of our control action for a predefined time horizon (see X line in Fig. 2.7.a). Second,
we solved the time-fractal optimal control problem to find the control signalu
fractal
(t)
and determined the level of the glucose for a predefined time horizon (see X line in
Fig. 2.7.a). As can be seen from Fig. 2.7.a the expected blood glucose level of the
conventional (integer order) controller is different from what the integer order controller
expects. In other words, due to the fact that integer order controller ignores the long
term memory nature of the blood glucose, the control signal solution it proposes does
not bring the blood glucose to the desired range properly.
As discussed in Section 2, the CPS process of sensing the blood glucose can be
affected by either measuring errors (e.g., physical pressure on sensors due to body pos-
ture leading to inaccurate glucose levels, parameter estimation error) or intrinsic vari-
ability (e.g., inhomogeneous insulin and meal absorption). Consequently, in what fol-
lows, we investigate the noise effect on controlling the blood glucose during the finite
control horizon. We assume that the noise can be statistically characterized by Gaussian
32
Figure 2.7: a) In-silico blood glucose level by applying integer and fractal optimal con-
trol signals. b) Histogram of blood glucose levels of the integer order controller due to
Gaussian perturbations. c) Histogram of blood glucose levels for the fractal controller
in the presence of Gaussian perturbations.
distribution with mean 0 and variance 5. Figures 2.7.b and 2.7.c illustrate the effect of
the Gaussian noise on both the non-fractal (integer order) and time dependent fractal
(fractional order) controller. As can be seen from these plots, due to the fact that non-
fractal controller cannot accurately capture the complexity of blood glucose dynamics,
the presence of noise can decrease the blood glucose at dangerous levels (i.e., below
80mg=dL) representing a life threatening situation. In contrast, the time dependent opti-
mal controller is able to keep the blood glucose dynamics within safe bounds. A similar
trend and better control performance for the fractal controller has also been observed for
noise sources drawn from other distributions (e.g., Gamma, Levy).
33
Chapter 3
Comprehensive dynamic analysis of
gene expression
3.1 Gene expression: definitions, modeling, and chal-
lenges
Protein synthesis is fundamental for biological systems to perform a variety of func-
tions. They control the organism’s shape or can function as enzymes catalyzing specific
metabolic pathways to regulate specific cellular processes. These cellular functions
include responding to stimuli, transporting molecules and catalyzing metabolic reac-
tions. In order to program cells for performing the desired functionality, one should
regulate the protein synthesizing process. The process of protein synthesis from the
activation of a specific gene is called gene expression [46, 47, 48, 4].
Gene expression (briefly shown in Figure 3.1.a) is the process in which the genetic
information of a cell causes a cell to generate a functional gene product and, finally,
perform specific cell functions [49]. In other words, it is the process by which genotype
information gives rise to phenotype (observable characteristics). It is a vital process,
which causes cellular differentiation, morphogenesis, and the versatility and adaptabil-
ity of any organism [50]. Controlling the production process of the desired gene expres-
sion product (e.g., a protein) refers to the regulation of the gene expression process. The
regulation of gene expression controls the amount and timing of production of target
34
proteins [51]. Hence, investigating the dynamics of gene expression enables to under-
stand the mechanisms driving biological organisms. This knowledge helps us from both
scientific and engineering perspectives. It can be exploited to detect an anomaly or dis-
ease or to engineer cells for performing specific tasks (e.g. drug delivery for cancer
treatment) as it is the target of synthetic biology.
The biophysical mechanism affecting the regulation process has been actively stud-
ied [52, 53, 54, 55, 56]. For instance, searching for the target gene by the transcription
factors is discussed in [52, 53] and the diffusion process of search for the target genes
is studied in [54, 55]. Also, the spatial distribution of gene products is reported in [56].
However, these prior studies were not concerned with the mathematical characterization
of the gene expression dynamics for several gene regulations in a network of genes. To
identify the main mathematical characteristics of gene expression dynamics, we inves-
tigate individual and cross-dependent gene expression time series. First, we investigate
the statistical properties of single (isolated) gene expression time series (shown in 3.1.c
and 3.1.d), and, then, we analyze the cross-correlation between pairs consisting of a
gene and a transcription factor in the gene regulatory network (Figure 3.1.b). In contrast
to the previous study [57] in which regulation of cell fate through genome-wide expres-
sion by temporal-spatial self-organization is considered, here, we mainly analyze the
temporal variability of individual genes. We then investigate the correlation of linked
transcription factor and genes. Moreover, we analyze the expression level of single cells
without considering population effect [58].
35
3.2 Mathematical definitions of cross-dependence frac-
tal dynamics
3.2.1 The Hurst exponent and Multifractal Detrended Fluctuation
Analysis (MFDFA)
We have used the MFDFA method for analysis of gene expression time series. This
method, which is the extension of detrended fluctuation analysis (DFA) to extract the
Hurst exponent [59], is introduced in [60] for analysis of multifractal properties of non-
stationary time series. Since we do not have the stationarity assumption for gene expres-
sion time series, MFDFA method is a suitable one for studying them. Scaling properties
and long-range dependency of time series can be obtained by the DFA method. How-
ever, for time series with multifractal properties and different scaling exponents on dif-
ferent scaling regimes or different time intervals, it is essential to exploit the multifractal
detrended fluctuation method (MFDFA) to reveal the multifractal property.
The MFDFA method consists of five steps to estimate the multifractal spectrum of
a nonstationary time series. Similar to the DFA method, the profile of the time series is
obtained first, which is determined by the integration of the difference of the time series
with its average value (^ x).
y(i) =
P
k
i=1
(x(i) ^ x)
Second, it divides the profile into non-overlapping segments (or scales (n= N/s),
where s is the scale). For each of these boxes, a least squared local trend is fitted.
Third, it calculates the local trend within each segment. For each of these boxes, a
least squared local trend is fitted. The value of the fitted time series obtained for boxes
of length (n) is denoted byy
n
.
36
Fourth, it computes the average of the fluctuation function over all segments to
obtain the qth order fluctuation function.
F (s)
q
=
P
Ns
k=1
((y(i)y
n
(i))
2
)
q
=N
s
Finally, the Generalized Hurst exponent is estimated by fitting a linear line to the
log-log plot of theF (s)
q
with respect to scale (s), according to the following equation:
F (s)
q
=s
H(q)
The Hurst exponent is the value of the Generalized Hurst exponent (H(q=2)), which
is a special case and is used usually when one is interested only in analyzing the long-
range dependency of a signal and not the multifractal chrematistics.
Finally, the multifractal spectrum of the multi-variable signal (;f()) is estimated
by the Legendre transform:
(q) =H(q) (q 1)
(q) =
d(q)
dq
f() =q(q)
3.2.2 Detrended Cross-Correlation Analysis (DCCA)
This method is designed to investigate the power law cross-correlation between two
time-series [61]. Similar to the detrended fluctuation analysis (discussed in the previous
sub-section), which computes the scaling behavior of the auto-correlation function, the
DCCA method computes the scaling behavior of the cross-correlation function between
two time-series and analyzes its scaling behavior.
DCCA method first computes the integrated profile of each time series:
37
y
1
(k) =
k
X
i=1
(x
1
(i) ^ x
1
)
y
2
(k) =
k
X
i=1
(x
2
(i) ^ x
2
)
Second, both the entire time series is divided into non-overlapping intervals.
Third, it computes the local trend in each interval for each time series
(y
1;n
(i)andy
2;n
(i)). Fourth, it calculates the covariance of the residual of profiles from
local trends. It calculates the detrended covariance (H(q)) by summing over all segments
of the nonstationary time series:
^
F (s)
q
=
Ns
X
k=1
((y
1
(i)y
1;n
(i)
2
(y
2
(i)y
2;n
(i)
2
)
q
=N
s
(3.1)
The cross-correlation exponent is estimated by fitting a linear line to the log-log
plot of the
^
F (s)
q
with respect to scale (s):
^
F (s)
q
=s
(3.2)
3.2.3 The Mandelbrot Binomial Cascade Model
This model is proposed by Mandelbrot [62] to better explain the erratic or fractal appear-
ance of a probability measure as an alternative for probability distribution. It starts with
a probability measure, which is self-similar:
([a;b]) =m
0
([2a; 2b]) +m
1
([2a 1; 2b 1]) (3.3)
38
Once the unit interval [0,1] is divided into two subintervals,m
0
mass is assigned to
the left subinterval and (m
1
= 1m
0
) is assigned to the right subinterval. Repeating
this step for each of the subintervals for n times will result in the Mandelbrot model
with n iterations. Mandelbrot has proved that the limit behavior of this model when n is
infinitely large (1) can be best illustrated by multifractal formalism. He formulated the
(;f()) spectrum on the basis of the parameters of the Mandelbrot Cascade model.
We have compared the observed multifractality spectrum in gene expression time series
to the closest one obtained by Mandelbrot cascade model.
3.2.4 Random Cascades on Wavelet dyadic trees
This model [63] is proposed to model multifractal objects. The notion of cascade here
refers to a self-similar process whose properties are defined multiplicatively in different
scales. In summary, in this model, the wavelet coefficients of a function are self-similar
at different scales. Two types of W-cascades are proposed: Log-Normal W-cascades
and Log-Poisson W-cascades. For the log-normal cascade (with and for parameters
of the normal random variable), the following equation holds for singularity spectrum:
F () =
ln2
2
2
( +
ln2
)
2
(3.4)
For the Log-Poisson W-cascades, the following equation holds for the mass expo-
nent:
(q) =
1
ln2
((1
q
)
q) 1 (3.5)
As can be seen, the second derivative of the mass exponent for Log-Poisson W-
cascade model has the following equation:
39
0(q) =
1
ln2
((1
q
))
ln2
00(q) =
ln2
(ln)
2
q
We have reported the similarity of the multifractal spectrum of cross-dependencies
in gene expression time series to the log-normal W-cascades model. Also, we have
reported the disagreement of the multifractal spectrum of cross-dependencies in gene
expression time series to log-Poisson W-cascade model due to its deviation from power-
law shape as shown in next chapter.
3.2.5 Entropy and Entropy of a Network
Shannon entropy [64] is a measure of the unpredictability of the state, or equivalently,
of its average information content. Shannon defined the entropy of a discrete random
variable X with possible values offx
1
;x
2
;:::;x
k
g and probability mass function P(X)
as:
H(X) =E[I(X)] =E[logP (X)] (3.6)
More explicitly, entropy can be written as:
H(X) =
n
X
i=1
P (x
i
)logP (x
i
) (3.7)
Entropy is a measure of the unpredictability of the state or its average containing
information. One example to illustrate is when there is no uncertainty and the random
variables take only one value in which the value of the entropy will be zero. As the
number of possibilities increases, the entropy increases as well.
40
We have used the notion of entropy in the context of networks. We consider the
weight of the links in the network as the random variable and we discuss entropy of
the weight of the weighted links. Given an undirected binary graph of gene regulatory
networks and time series of genes and transcription factors in the network (which are
the nodes in the gene regulatory network), we construct a weighted network (shown in
Figure 3.2). In the constructed weighted network, the weight of each link is the cross-
correlation exponent of the time series of two time-series linked together in the gene
regulatory network. Then, in the new constructed weighted network, we consider the
distribution of the weights of the links and entropy of them as a measure of the entropy
of the network. Also, this weighted network can be used for other algorithms measuring
the entropy of complex networks proposed in [65, 66]. Figure 3.2 illustrates this method
by showing how the weights are assigned to each link. In this Figure,
1;2
;:::;
1;9
are the
cross-corrections of the time series of the transcription factor and genes which are linked
together in the gene regulatory network in the left part of the figure. Hence, this shows
how knowing the existing interactions in the network and having the time series of each
nodes dynamics can lead us to know cross-correlation exponents and then assigning the
concept of entropy to the network dynamics.
3.3 Mathematical definitions of cross-dependence frac-
tal dynamics
3.3.1 Gene expression dynamics exhibits long-range dependency
and multifractal properties
We investigate the statistical properties of gene expression data and compute the Hurst
exponents of gene expression time series in Saccharomyces cerevisiae (S. cerevisiae)
41
and Escherichia coli (E. coli) . Figure 3.3.a and Figure 3.3.b show the log-log plot of
the fluctuation function as a function of the scale for the time series of a transcription
factor (ynel) in S. cerevisiae and E. coli, respectively. In these plots, the slope of the
curve represents the Hurst exponent. We observe that 95% and 98% of the time series
of genes from the S. cerevisiae and the E. coli gene expression networks, respectively,
exhibit a long-range dependency property. More precisely, their Hurst exponent was
greater than 0.5. To demonstrate this important property, Figure 2c and Figure 2e show
the histogram of the Hurst exponent of gene expression time in S. cerevisiae and E.
coli, respectively. Generally speaking, a Hurst exponent that exceeds the 0.5 threshold
value denotes a persistent (positively correlated) behavior in the sense that a high value
is likely to be followed by another high value with nonzero probability. In addition,
because the Hurst exponent for most of the genes is significantly higher than 0.5, the
gene and transcription factor time series cannot be regarded as a random process and
modeled through Markovian formalism [67]. This mathematical characteristic provides
a clue as to how to construct stochastic models for gene expression processes, but this
is left for future work. We observe the same property in the time series of transcription
factors inS.cerevisiae andE.coli. More precisely, 97% of the transcription factors inS.
cerevisiae andE.coli possess the long-range dependence property. The histogram of the
Hurst exponent of transcription factors in S. cerevisiae and E. coli are shown in Figure
3.3.d and Figure 3.3.f, respectively.
Employing fractal analysis is also helpful to gain insight into other interesting prop-
erties. Here, we see a bimodal characteristic in the Hurst exponent distribution, shown in
Figure 3.3.c, Figure 3.3.d, Figure 3.3.e, and Figure 3.3.f. This feature is especially vis-
ible in Figure 3.3.c and Figure 3.3.e, where the histogram of genes in S. cerevisiae and
E. coli is presented. This may be explained by a possible bimodal diffusion potential,
as in [68]. Since gene expression includes a diffusion process with multiple diffusion
42
potentials (inside and outside the nucleus), this bimodality can be explained by non-
equilibrium Brownian motion with multiple potential profiles. However, further exper-
imental studies are required to elucidate the nature and implications of these bimodal
statistics.
By employing the multifractal detrended fluctuation analysis (MFDFA) [60] (see
Methods section) to examine the multifractal characteristics of gene expression time
series, we observe that both genes and transcription factors have pronounced multi-
fractal properties. For monofractal behavior, the generalized Hurst exponent displays
a linear dependency with the order q of the cross moments. Instead, if the general-
ized Hurst exponent exhibits a nonlinear dependency (such as the S-shape displayed
in Figure 3.4.b) as a function of the order q of the cross-moments, then the stochastic
interdependence is considered to possess multifractal characteristics.
To provide a more in-depth report, we use the bootstrapping technique [69] to inves-
tigate the existence of the long-range dependence property, considering the limitations
related to the length of the experimental time series. For every gene expression time
series, we have sampled 10 randomized subintervals of the gene expression time series,
each containing 90% of the ordered piece of the original time series. Then, we calcu-
late the Hurst exponents for all the randomized versions. The difference between the
percentage of the long-range dependency for the gene expression time series and the
randomized versions was approximately 0.006 forE.coli and 0.0001S.cerevisiae.
43
3.3.2 The distribution of cross-correlation exponents of pairs of
genes and transcription factors of gene regulatory networks
has a wide range of variation
Although we observe the fractal and long-range cross-correlation in linked pairs of
genes and transcription factors in the gene regulatory networks, the cross-correlation
exponents were not the same in all the links. We have shown the distribution of the
cross-correlation exponents for pairs of genes and transcription factors in the S. cere-
visiae and E. coli gene regulatory networks in Figure 3c and Figure 3d, respectively.
Inspired by Shannon entropy [64], we use this histogram for measuring the entropy, and
hence, the information content of a gene regulatory network across different cell types
for quantitative analysis and specification of gene regulatory networks. The computed
Shannon sample entropy for emph S. cerevisiae and E. coli was 4.18 and 5.29, respec-
tively. Consequently, we conclude that the gene expression network ofE.coli has more
complex dynamics than that of S. cerevisiae. Also, considering a static gene regulatory
network and having traces of gene expression time series for a cell at different times, we
can compute the cross-correlation exponents for the links at a different time. This can
be useful to compare statistical properties and complexity of dynamics. Similarly, by
having different time series of gene expression dynamics, we can compare normal vs.
disease affected (for example cancer typed) cells.
3.3.3 Multifractal characteristics of interactions within the gene
regulatory network can be modeled by random cascades on
wavelet dyadic trees
We analyzed the multifractal property of the cross-correlation of pairs of genes and
transcription factors in a gene regulatory network. We investigated whether the observed
44
multifractality can be explained by the known analytical cascade models including the
Mandelbrot bimodal cascade model [62] (see Methods section) and the random cascades
on wavelet dyadic trees [63]. We observe deviations of the empirical spectrums from
the Mandelbrot model and an approximate agreement to the random cascades on wavelet
dyadic trees model.
Based on the range of the Holder exponent values in the multifractal spectrum, we
observe that only 0.04 of the links in a gene regulatory network of S. cerevisiae and
none of the links in the network of E. coli can be modeled by the Mandelbrot cascade
model for multifractal spectrums (see Methods section). We observe that even for the
few links that we could find a closest Mandelbrot model spectrum, the deviation from
the Mandelbrot model and the data we had for gene regulatory network was significant.
We show two such samples in Figure 4. Figure 4a shows several multifractal spectrums
for the links in E. coli gene networks. Note that the peak of the multifractal spectrum
for these spectrums was lower than the value 1, which does not fit with the Mandelbrot
Binomial Cascade Model [62]. Figure 4b shows several multifractal spectrums for the
links inS.cerevisiae. Figure 4c shows the closest Mandelbrot Model we could fit for the
links in theS.cerevisiae gene regulatory model. There is a significant deviation between
the Mandelbrot model and the spectrum from gene regulatory network data.
We also investigated the agreement between the observed multifractality of the
cross-dependencies in the gene regulatory network of S. cerevisiae and E. Coli, respec-
tively, and a few well-known multifractal models such as the random cascades on
wavelet dyadic trees [63]. We investigated whether the log-normal W-cascade model
can be fitted to the cross-dependencies (links) in the two above-mentioned gene reg-
ulatory networks. We extracted the parameter of the estimated log-normal W-cascade
model based on the peak of the empirical spectrum and the variation of the singularity
45
spectrum (see the Methods section). We observed very similar spectrums for a signif-
icant number of links. In this study, we used two threshold values of 70% and 75%
for the ratio of mentioned areas. For the gene expression cross-dependencies in E. coli,
we observed a 74% and 38% agreement between the postulated model and the empiri-
cally estimated multifractal spectrums when considering overlapping area ratio thresh-
olds of 70% and 75%, respectively. For the gene expression cross-dependencies in S.
cerevisiae, we observed a 59% and 31% agreement between the postulated multifractal
model and the empirically estimated multifractal spectrums when considering overlap-
ping area thresholds of 70% and 75%, respectively. Figure 5a and 5b show a best fitting
scenario between the postulated multifractal model and an empirically estimated spec-
trum for a cross-dependence in the gene regulatory network ofE.Coli andS.cerevisiae,
respectively.
We also investigated the agreement between the empirical multifractal spectrums
and the log-Poisson W-cascade model [63]. We observed that the empirical multifractal
spectrums could not be described by this cascade model since the second derivative of
the mass exponent should follow a power law (see Methods section) while the empirical
data has a significant deviation from a power law trend (see Figure 3.7).
3.4 Implications of cross-dependence fractal behavior
in gene expression
A genome expression vector is the most informative descriptor of a cell state, as the
functional state of an organism is determined largely by the pattern of expression of its
genes [5]. Gene expression is the process in which information from a gene is used to
synthesize a functional gene product. It is the process in which the information flows
46
within a complex biological system. As the search for patterns in nature and their inter-
pretation is one of the main purposes of science, unveiling the DNA patterns in those
sequences has become an exciting challenge to the present generation of biologists, sta-
tistical physicists, and information scientists. Toward this end, many researchers have
studied the statistical properties of coding and non-coding segments of DNA sequences.
They have reported interesting results showing fractal nature of coding DNA regions
[46, 47, 48, 4]. However, these studies fail to address the dynamical properties of the
biological systems. Since biological systems are dynamic, their study requires monitor-
ing their activity at multiple time points.
To investigate the causal relations in gene expression, numerous biophysical mecha-
nisms affecting the regulation process were studied in [52, 53, 54, 55, 56]. It is demon-
strated by several simulations that rapid and reliable gene regulation requires that the
transcription factor (TF) be close to their target site on DNA [52]. In [53], the authors
use an explicit model for numerical analysis. The authors report that the observed varia-
tions in regulation efficiency are linked to the magnitude of the variation of TF concen-
tration peaks as a function of binding site distance from the source. In [54], the authors
have presented a semi-analytical model for the in vivo target search of the TFs within a
diffusion framework. They have shown that alternating between three-dimensional bulk
diffusion and one-dimensional sliding along the DNA contour can provide a quantitative
approach to gene regulation in living bacteria cells. Their proposed model agrees with
experimental findings regarding the mean search time of lac repressor in a living E. Coli.
In [55], the authors have reported their observation of kinetics of the gene expression
process at the single-molecule level in living cells by labeling with a fluorescent pro-
tein, which agrees with 1D diffusion along DNA segments and 3D diffusion. In [56],
the authors test the expectation of previous theoretical models using high-throughput
47
single-molecule microscopy to determine the average spatial distribution of lac repres-
sor. Their finding shows inconsistency between expectations and experimental findings.
They show that the gene products distribution is spatially inhomogeneous and depen-
dent on the location of the repressor gene in bacteria and eukaryotes. However, they do
not consider the gene expression dynamics from a network of genes perspective and do
not account for cross-correlations and multiscale phenomena.
To study the dynamic nature of gene expression processes, researchers must moni-
tor activity levels of genes and transcription factors at multiple time points. The most
informative source of information regarding gene expression activity is the gene expres-
sion time-series. With advances in gene array technology, the level of gene expression
of thousands of genes (by providing the concentration level of gene expression prod-
ucts) can be measured simultaneously. By accessing a high- throughput data collection,
a wide range of insights, such as characterizing the functions of specific genes, the
relationships among these genes, and their regulation and coordination can be gained.
These insights can also be used to understand the gene regulatory network as a complex
network. There are many studies which try to infer the underlying gene regulatory net-
work from empirical time series [5]. However, little is known about the mathematical
characteristics of the gene expression dynamics from a complex systems perspective.
In this paper, we investigated the scaling properties of gene expression dynam-
ics. Unlike previous work that demonstrated the fractal properties in DNA sequences
[70, 71, 72, 5], we investigate the dynamics of cross-dependencies between genes and
transcription factors within the gene regulatory networks. We show that the gene expres-
sion time series (which is the concentration of gene expression products in the process
of gene expression) have fractal and long-range dependence properties in E. coli and
S. cerevisiae. We also investigate the cross-correlation of gene-transcription factors,
48
which are linked together in gene expression networks. We report the fractal and long-
range cross-dependency of linked genes and transcription factors of gene expression
networks in E. coli and S. cerevisiae. We also show that the multifractal nature of these
cross-correlations cannot be modeled through a Mandelbrot binomial cascade model. In
contrast, we found very good agreement between the empirical multifractal spectrum of
the cross-dependencies in the gene regulatory networks and the log-normal W-cascade
model. We suggest investigating more cascade models on empirical data [73, 74, 75]
as future work. In summary, there is a need for more advanced theoretical models that
can capture the multifractality observed in this critical biological process. One possible
method for modeling gene expression dynamics can exploit the multifractal Fokker-
Planck formalism, as discussed in [76].
We also propose using the distribution of cross-correlation exponent of the links in
gene regulatory network as a measure of the complexity for the regulatory networks.
Having this complexity measure enables a quantitative descriptor for different cell types
or to differentiate different cell fates when the system undergoes transitions. We report
the distribution of cross-correlation exponent of links in regulatory networks of E. coli
and S. cerevisiae as case studies. We suggest investigating this property on a wider range
of biological systems when enough data sets are available. We also propose using this
property as a network property in general. We propose using the distribution of cross-
correlation exponents of gene-transcription factor links in a complex network to measure
the complexity of the interactions in the network. Also, the computed cross-correlation
exponents of a network can be used by other algorithms, such as those proposed in
[65, 66, 77] for computing the entropy of a network. This has a variety of applications
in distinguishing the different status of cells (e.g. healthy vs. disease affected states).
This can reveal insightful results in many complex networks either in biology, social,
financial and many other interesting examples of complex networks.
49
Our findings explain the inherent variability in gene expression processes, even
among isogenic cells situated in an identical environment. Because of the long-range
cross-dependency of a gene and its linked transcription factor, the current concentration
level of a gene depends both on current and previous values of its own and its linked
transcription factor. As explained in previous studies, this leads to phenotypic diver-
sity, which can be helpful for surviving in an uncertain and fluctuating environment
[78]. Also, the endogenous cellular mechanism through positive and negative feed-
back controls variability in gene expression to prevent disruption from normal devel-
opment. Hence, unlike the usual assumption about noise as a nuisance, variability in
gene expression makes the population of cells more robust against environmental fluc-
tuations. Interestingly, there are other examples in nature in which the presence of noise
makes the system smarter. For example, in [79], they have shown how the presence of
noise in a network of spiking neurons in the human brain enables probabilistic reasoning
and creative problem-solving.
This study is the first to demonstrate the long-range dependency of gene expression
dynamics. In contrast to previous studies [80], which have shown the long-range depen-
dency for the structure of DNA, we investigate the dynamics of gene expression time
series. Previous studies show that coding regions of DNA structure, which store the
biological information for the gene expression process, possess the long-range depen-
dency property [80]. In contrast, our results report the same property in gene expression
time series. Of note, these dynamics stem from the transformation of information from
the structure to dynamics by producing gene products. This is an insightful empirical
result that can trigger more studies on other examples from nature, as well as analytical
and mathematical investigations. For example, investigating other processes that follow
a rule from a static structure to generate dynamical products and process (such as the
50
central dogma of biology) can be interesting. Lastly, mathematical and analytical inves-
tigation of the relation between structure and dynamics of processes are also fundamen-
tal in theory. It would be revealing to investigate how long-range dependency (and/or
fractal/multifractal properties) evolves from structure to dynamics (and vice versa) in
processes. Answering to the question of how long-range dependency transfers between
structure and dynamics and how the degree of fractality/multifractality of structure and
dynamics are like each other would have a huge impact on predicting the behavior of
complex systems.
51
Figure 3.1: (a). Gene expression process is initiated with the triggering of a gene in
the DNA strand and continues with the generation of RNA, mRNA and, finally, the
protein product. (b). Part of the gene regulatory network of E. coli reported by [4],
which consists of genes and transcription factors. The diamonds represent transcription
factors and the circles represent genes. (c). Time series of transcription factors of E. coli
obtained from [5] and explained in the methods section. (d). Time series of genes of E.
coli obtained from [5] and explained in the methods section.
52
Figure 3.2: Extracting a weighted graph from cross-correlation exponents of the linked
nodes in a gene regulatory network
Figure 3.3: (a). Scaling of fluctuation function of a gene time series in S. cerevisiae,
(b). Scaling of fluctuation function of ynel gene time series in E. coli, (c). Histogram of
cross-correlation exponent of genes in S. cerevisiae, (d). Histogram of cross-correlation
exponent of transcription factors in S. cerevisiae, (e). Histogram of cross-correlation
exponent of genes in E. coli, (f). Histogram of cross-correlation exponent of transcrip-
tion factors in E. coli.
53
Figure 3.4: (a). Scaling of detrended covariance function in a gene-transcription factor
link in E. coli (ihfB to ompR) and S. cerevisiae (YLR256W to YKL020C) regulatory
network, (b). Generalized cross-correlation exponent of a gene-transcription factor link
in E. coli and S. cerevisiae, (c). Histogram of cross-correlation exponent of genes in
S. cerevisiae, (d). Histogram of cross-correlation exponents of gene regulatory network
links.
54
Figure 3.5: (a). Multifractal spectrum of several randomly picked links of gene regu-
latory network of E. coli, (b). The multifractal spectrum of several randomly picked
links of gene regulatory network of S. cerevisiae, (c). The multifractal spectrum of a
link between the gene regulatory network of S. cerevisiae and the best spectrum from
the Mandelbrot binomial model.
Figure 3.6: (a). The multifractal spectrum of a cross-dependency (link) between the
gene regulatory network of E. coli (soxR and soxS genes) that can best be fitted by
the log-normal cascade multifractal model, (b). The multifractal spectrum of a cross-
dependency (link) between the gene regulatory network of S. cerevisiae (YKL032C and
YKL043W genes) that can best be fitted by the log-normal cascade multifractal model.
55
Figure 3.7: (a). The second derivative of the mass exponent for cross-dependencies in
E. coli. (b). The second derivative of the mass exponent for cross-dependencies in S.
cerevisiae.
56
Chapter 4
Mutli-fractal multi variable detrended
fluctuation analysis
4.1 Medical Complex (biological) systems: definitions
and challenges
Elucidating the dynamic nature of complex systems has important consequences not
only on deciphering the topological interactions among various system components but
also for predicting their behavior. To date, most of the mathematical modeling research
for complex systems have ignored the spatio-temporal (multi-scale) correlations and
modeled the complex systems through memoryless (integer order differential or differ-
ence) equations [81, 82]. Simply speaking, the current mathematical models of complex
systems neglect either the observed fractality in the dynamics of individual components
of the complex system or the higher-order cross-correlations between various compo-
nents of the complex system [83, 84, 85, 86]. To remedy the disconnect between the
current mathematical modeling [87, 7, 88, 89] and the intrinsic mathematical character-
istics of complex systems, in this paper, we introduce a multi-dimensional (-variable)
multi-fractal detrended fluctuation (mMFDFA) analysis method that allows to quantify
the cross-correlation complexity and memory existing between a group of interacting
components of a complex system. More precisely, our proposed method analyzes an
57
approximation of the higher-order statistics of multi-dimensional signals and quanti-
fies the power law scaling exhibited by dynamical fluctuations of coupled processes
(time series encompassing the evolution of the complex system components). To illus-
trate the depth and benefits of this mathematical method, we investigate the multi-fractal
spatio-temporal characteristics of a human brain. We analyze the spatio-temporal depen-
dencies among groups of neurons within the human brain during various sleep or task
performing periods. We demonstrate that the coupled brain dynamics exhibits multi-
fractal cross-correlations with specific patterns in terms of shape and wideness of the
multi-fractal spectrum. In addition, the proposed mMFDFA helps us quantify the multi-
fractal cross-correlations existing between human brain signals and physiological pro-
cesses. We observe the existence of different levels of multi-fractality across various
brain activities and sleep stages. The dynamics of brain activity is of interest for sev-
eral communities [90, 91]. Several research studies showed that individual brain signals
exhibit a multi-fractal behavior over time [92, 93]. In contrast to this body of work,
we analyze holistically (as an interdependent complex system) both the dynamics of
brain activity and the coupled dynamics between the brain and physiological processes.
Our research shows not only that the brain activity and physiological processes dis-
play a complex spatio-temporal multi-fractal behavior, but also that the degree of multi-
fractality is richer (wider) for more complicated functionalities in brain activities.
58
4.2 Mathematical description of multi-fractal multi-
variable detrended fluctuation analysis
4.2.1 Description of the method
The modified multi-variable multifractal detrended fluctuation analysis builds on clas-
sical single variable MFDFA34 with an extension in the definition of the fluctuation
function:
Step 1: Given M variables x
1
(i);x
2
(i);:::;x
M
(i) each having N number of sam-
ples, this method starts with constructing M time seriesy
1
(i);y
2
(i);:::;y
M
(i) with the
following definition:
y
j
(i) =
k
X
i=1
(x
j
(i) ^ x
j
) (4.1)
Where ^ x
j
=
1
N
P
N
i=1
x
j
(i) is the average value ofx
j
(i). Please note that we do not
assume anything about these variables having independent distribution for now and we
leave the discussion in numerical simulation section. The next step is to measure the
scaling characteristics of the signals:
Step 2: To measure the scaling behavior of the fluctuation of time series, each time
series is divided into boxes of equal length n (n =N=s, where s is the scale). For each
of these boxes, a least squared local trend is fitted. The value of the fitted time series
obtained for boxes of length (n) is denoted by ^ y
j
. For a given box size n, the average of
the qth order fluctuation of each time series (y
j
) with respect to the trend ( ^ y
j
) is:
59
F
2
j
(v;s) =
1
s
s
X
i=1
f ^ y
j
(i +sv)y
j
(i +sv)g
2
F
q
j
(s) =f
1
N
N=s
X
v=1
F
2
(v;s)
q
=2g
1=q
To consider the fluctuation of the overall multi-variable time series, consider the
product of the individual fluctuation functions:
d
F (s)
q
=
M
v
u
u
t
M
Y
j=1
F
j
(s)
q
(4.2)
Step 3: To calculate the scaling behavior of the overall fluctuation function, estimate
the Generalized Hurst exponent of the multi-variable signal by fitting a linear line to the
log-log plot of the
^
F (s)
q
with respect to scale (s):
d
F (s)
q
=s
H(q)
(4.3)
Step 4: Estimate the scaling exponent ((q)):
(q) =H(q) (q 1) (4.4)
Step 5: To characterize the multi-fractal formalism by the singularity spectrum
(f()), we use the following equations that relatef() to(q) via Legendre transform:
(q) =
d(q)
dq
f() =q(q)
60
Where is the singularity strength or Holder exponent andf() denotes the dimen-
sion of the subset of the series that is characterized by.
4.2.2 Computation complexity of mMFDFA method
Our proposed method has a computation complexity that scales linearly with the number
of variables of the multi-variable data (M) and the number of samples (N) of the data.
Since the main computation of the algorithm is the linear regression, the overall compu-
tational complexity of the mMFDFA method isO(NM). There is a constant number
of scales over which detrending is computed and that does not affect the computation
complexity.
4.2.3 The correctness of the proposed Algorithm
The main idea behind the proposed algorithm is on the basis of the fact that each of the
individual time series has an individual power-law fluctuation function:
F (s)
(q;i)
s
H
i
(q)
(4.5)
Hence, the product of their fluctuation function also has the power-law scaling
behavior:
d
F (s)
q
=
M
v
u
u
t
M
Y
j=1
F
j
(s)
q
=
M
p
s
H
1
(q)
:::s
H
M
(q)
=
M
p
s
H
1
(q)+:::+H
M
(q)
=s
(
P
M
i=1
H
i
(q)
M
)
Also, the generalized Hurst exponent of the multi-variable time series is the average
of the individual Hurst exponents of individual time series:
61
d
F (s)
q
s
(
P
M
i=1
H
i
(q)
M
)
(4.6)
The above algorithm is consistent with previously proposed algorithms for comput-
ing detrended cross-correlation exponent [94].
There are other algorithms [6, 95] that try to address the Generalized Hurst exponent
of multi-variable time series. However, instead of using the product of the individual
fluctuation function (shown in equation (3)), they use vector norm. We have argued
analytically in section 2.2 that by using our method and the assumption of power-law
relation for individual time series, the overall fluctuation function
d
F (s)
q
, will also have
power-law relation with respect to scale. In contrast, there is no guarantee that the sum of
two or more power-law functions (s
(h
1
(q))
+s
(h
2
(q))
++s
(h
M
(q))
as also power-law relation
with respect to scale. Although they may use the best fitting in the log-log plot, there
is no strong analytical proof that the norm obeys the power-law relation with respect to
scale. In the next section we provide numerical simulation results and in section IV we
provide results of the application of this method on human brain activities.
4.3 Verification of multi-fractal multi-variable
detrended fluctuation analysis on synthetic case
studies
We have validated the proposed mMFDFA on the multi-variable data using different
situations: On correlated time series. In these conditions, we make N = 216. We have
also considered the application of the mMFDFA method on independent time series
ranging from 2 to 100 variables.
62
4.3.1 Application of the algorithm on independent time-series
If a set of multivariable time series (X = x
1
(i);x
2
(i);;x
M
(i)) consists of individual
multifractal time series (x
1
(i);x
2
(i);;x
M
(i)), then the individual fluctuation functions
for each variable have a power-law relation with scale s:
F
(q;i)
(s)s
H
i
(q)
i
(q) =H
i
(q) (q 1)
In which the individual mass exponents,
i
(q), have a linear relation with respect to
q. By looking at equation (1) and (2), we can see that the multi-variable Generalized
Hurst exponent and mass exponent ((q)) will also have a linear relation with respect
to q. In conclusion, if the individual time series of multi-variable data are mono-fractal,
then by applying the proposed algorithm, the multi-variable multi-fractal analysis leads
to mono-fractal spectrum.
Figure 4.1, shows the scaling behavior of the fluctuation function for a bi-variable
time series consisting of two non-identical independent time series in Figure 1.a and
100-variables in Figure 1.b. As one can observe from Figure 4.1.a, when consider-
ing the bi-variable case of two time-series which are individually characterized by a
Hurst exponent of 0.7 (time series 1 in blue) and 0.9 (time series 2 in red) respectively,
the resulting Hurst exponent for bi-variable case is 0.8 (shown in green) which is the
average of the Hurst exponent of individual time series. We have also shown the multi-
variable multi-fractal spectrum of time series consisting of 10, 20, 30, 50 and 100 time
series in Figure 4.1.b. As can be seen, we see a shrinking trend in the multi-fractal
spectrum as we increase the number of variables taken into consideration. This comes
from the averaging nature of multi-variable Generalized Hurst exponent. Since these
63
Figure 4.1: Simulation Results on independent time series:(a). Detrended fluctuation
function of the independent bi-variable signal along with detrended fluctuation analysis
of individual time series (b). The multifractal spectrum of independent time series (we
have shown for bi-variable, 10, 40, 60 and 100 independent variables)
time series are independent, they do not get their maximum or minimum value of Hurst
exponents at the same time. On the contrary, for correlated time series, since the individ-
ual time series get maximum (or minimum) of their Hurst exponents at the same time,
the multi-variable General Hurst exponent will get the maximum (or minimum) range of
the average of individual time series. We have provided these results in the subsection.
4.3.2 Application of the algorithm on correlated time-series
If we apply the mMFDFA on two correlated time-series (x andy) which are constructed
using surrogate time seriesz (x andz are independently generated):
y(n) =x(n) +
p
1
2
z(n) (4.7)
Since each of thex andy time series, are having a Hurst exponent (H
1
(q) andH
2
(q)):
64
F
(q;1)
(s) =s
H
1
(q)
;F
(q;2)
(s) =s
H
2
(q)
(4.8)
from equation (3) of Step 2 of the mMFDFA algorithm we have:
d
F (s)
q
=
p
s
H
1
(q)
s
H
2
(q)
=s
((H
1
(q)+H
2
(q))=2)
(4.9)
Thus, the combined Generalized Hurst exponent and the multi-fractal spectrum of
two variables will be average of two of them.
If we continue increasing number of correlated variables by using equation (15), we
get the same result for Generalized Hurst exponent of the multi-variable signal:
d
F (s)
q
=
p
s
H
1
(q)
s
H
2
(q)
:::s
H
M
(q)
=s
P
i=1
M
H
i
(q)
(4.10)
In Figure 4.2.a, we have shown the log-log plot of the detrended fluctuation with
respect to scale for individual and combined time series for two correlated time series.
In this picture, time series x has the Hurst exponent of (H1=0.80) and time series y has
the Hurst exponent of (H2=0.72) and the Generalized Hurst exponent of the bi-variable
time series is 0.76 as we expected.
In Figure 4.1.b, we have shown the multifractal spectrum of correlated time series.
We have also shown the multi-fractal spectrum for 100 correlated time series we gen-
erated by ARFIMA process [96]. We have shown the results for the bi-variable case
as well as 10, 40, 60 and 100 correlated time series. As can be seen, the multifrac-
tal spectrum of the multi-variable time series includes the spectrums of the containing
multi-variable time series. Hence, the multi-variable spectrum of 100 correlated time
series (shown in red) contains the spectrum of 60 correlated time series (note that these
60 individual time series s are included in 100 correlated time series). As we increase
65
Figure 4.2: Simulation Results on correlated time series: (a). Detrended fluctuation
function of the correlated bi-variable signal along with detrended fluctuation analysis of
individual time series (b). The multifractal spectrum of correlated time series (we have
shown for bi-variable, 10, 40, 60 and 100 correlated variables)
the number of variables, we see the inclusion nature of multi-variable multi-fractal spec-
trum more obvious. This behavior comes from the averaging nature of Generalized
Hurst exponent and the fact that individual correlated time series get their minimum and
maximum of their Hurst exponent at the same time. Hence, the average of them gets a
bigger range compared to when these individual Hurst exponents are independent.
For the purpose of comparison, we have compared the results of our method to the
method proposed before [6, 95] in Figure 4.3. Figure 4.3.a shows the multivariable
multi-fractal spectrum obtained by our proposed method for up to 100 correlated time
series generated by ARFIMA [96] process and Figure 4.3.b shows the multivariable mul-
tifractal spectrum obtained by the previous method proposed in [6, 95]. As can be seen,
in contrast to our method, the previous method does not capture the inclusive nature of
multi-fractal spectrum for the correlated multi-variable time series. Therefore, on the
basis of the analyses in Section III, we can generally draw to the following conclusions:
(1) for the multivariate time series, the scaling exponent H obtained by MMFDFA is
66
Figure 4.3: The multifractal spectrum of correlated time series by different meth-
ods: a). The multifractal spectrum of correlated time series (we have shown for a
bi-variable case, 10, 40, 60 and 100 independent variables) obtained by our proposed
MMFDFA, (b). Multifractal spectrum of correlated time series obtained by the previous
method proposed in [6, 7]
basically identical with the average of the individuals if all the initial data channels are
generated from the same system whether they are mutually independent or correlated;
(2) For correlated multivariable time series, the multifractal spectrum obtained by our
proposed method includes the multifractal spectrum of individual time series. We have
shown that this property is not true for previously proposed method [6, 95] for mul-
tivariable time series and they could not distinguish between correlated multi-variable
time series and independent ones; (3) For independent multivariable time series, the
multivariable spectrum, does not necessarily include the individual time series.
67
4.4 Verification of multi-fractal multi-variable
detrended fluctuation analysis on brain datasets
Understanding the dynamics and complexity of the human brain is fundamental not
only for developing new medical therapies but also for developing new potentially more
powerful and efficient computing paradigms and deciphering the emergence of cog-
nitive (thinking) activities. Consequently, there have been numerous graph-theoretic-
based approaches to characterize the brain structure and measuring its properties such
as the small-world and self-similarity [97, 98]. Several works also studied the behav-
ior and dynamics of the brain signals in isolation as well as the functional interdepen-
dences between brain and body [99, 100]. Sleep plays an important role in resting the
whole body and dreaming is essential for the development of mental abilities and cre-
ativity. Several research efforts investigated and characterized the difference between
the human brain dynamics in sleep (considering different sleep stages) and wake states
[101, 102, 103]. Most of these prior research efforts ignore the complexity of the body-
brain interaction and fail to address the collective dynamics of brain-body activities.
Using the mMFDFA method, we measure the multi-fractal spectrum of brain-body inter-
action. To discriminate between different sleep stages and to compare the brain activity
of different subjects, we quantify the complexity by the strength of a multi-fractal spec-
trum. More precisely, we calculate the multi-variable multi-fractal spectrum and use
the wideness of the multi-variable multi-fractal spectrum as a measure of complexity.
Simply speaking, we say that a system is more complex than another one if its multi-
variable multi-fractal spectrum is wider. To the best of our knowledge, this is the first
attempt to characterize the complexity of brain-body interaction system and collective
dynamics of brain activities. We use a dataset [104] consists of EEG signal measure-
ments from several patients. We perform the mMFDFA method on four EEG signals
68
Figure 4.4: (a.) The multifractal spectrum of EEG time series for a healthy person (n1)
in S3 deep sleep state (b). The multifractal spectrum of EEG time series for a healthy
person (n1) in the REM sleep state (c). The multifractal spectrum of EEG time series
for a healthy person (n2) in S3 deep sleep state (d). The multifractal spectrum of EEG
time series for a healthy person (n2) in a REM sleep state
(Fp2-F4, F4-C4, C4-P4, P4-O2). We show that by applying the mMFDFA method we
can better distinguish among different states of the brain when performing various sleep
stages.
69
4.4.1 Even a sleeping brain exhibits a complex coupled multi-
fractal dynamics
We employ the mMFDFA method on a coupled sleeping brain-body dynamic repre-
sented as a 4-channel signal consisting of 4 EEG signals, namely Fp2-F4, F4-C4, C4-P4,
and P4-O2 for each patient. To contrast the multi-fractal behavior between individual
EEG and the coupled brain-body dynamics, Figure 4.4 shows the Multi-variable Multi-
fractal spectrum for the following cases: Figure 4.4.a and 4.4.b show the Multi-variable
multifractal spectrum of individual and combined EEG signals for a healthy individual
for S3 deep sleep stage (in Figure 4.4.a) and Rapid Eye Movement (REM) (in Figure
4.4.b). Figure 4.4.c and 4.4.d show the Multi-variable multifractal spectrum of individ-
ual and combined EEG signals for another healthy individual for S3 deep sleep stage
(in Figure 4.4.a) and Rapid Eye Movement (REM) (in Figure 4.4.b). As can be seen,
even sleeping brain has multi-fractal dynamics both for individual EEG signals as well
as combined multi-variable signals.
4.4.2 The multi-variable multi-fractal analysis reveals the distinc-
tion of dynamics of deep sleep stages
By employing our mMFDFA method on coupled sleeping brain-body dynamics, we
realized that multi-variable multifractal spectrum for Rapid Eye movement (REM) stage
has lower complexity compared to S3 stage of the sleep. The multi-variable multi-fractal
spectrum of REM stage was significantly narrower than that of S3 for all 10 health sub-
ject studies in [104]. As can be seen from Figure 4.4, this shortening in the multi-fractal
spectrum is not visible for individual EEG signals. However, the combined 4-EEG sig-
nal of brain signal had different multi-fractal spectrum wideness. This can be explained
by our findings regarding correlated signals. In other words, there is more correlation
70
among individual EEG signals when in S3 deep sleep stage compare to REM stage. In
the REM stage, EEG signals have less correlation and there is shrinking behavior in the
coupled EEG signals.
4.5 Verification of multi-fractal multi-variable
detrended fluctuation analysis on river network
datasets
4.5.1 Water flow dynamics across large river networks exhibits
multi-fractal characteristics.
In contrast to previous work that investigated the river dynamics for individual river flow
data, we applied the mMFDFA method concomitantly to multiple water flow time series
across the river basin. We used water flow data from the upper Mississippi (see Figure
4.5.a) and the Sacramento river basins (see Figure 4.5.b). Figure 4.5.c and 4.5.d show the
river flow dynamics at different sites across these river basins during October 2014 (i.e.,
the Mississippi river flow levels measured at Royalton, St. Cloud and Brooklyn Park
in Minnesota and the Sacramento river flow levels measured at Keswick, Red Bluff and
Grimes in California). Figure 4.5.e and 4.5.f show the generalized Hurst exponent as a
function of the q-th order of the statistical moments (of fluctuations) for both individual
river channels and the multi-variable integration of river channels. One can notice from
Figure 4.5.c that while the Mississippi river flow levels exhibit similar trends their indi-
vidual and coupled multi-fractal characteristics are different. Moreover, we also notice
that the wideness of the generalized Hurst exponent is higher when analyzing each indi-
vidual river flow dynamics (i.e., the red, blue and green curves) than when analyzing
71
the combined multi-variable dynamics (black curve). A similar shrinking trend in the
degree of multi-fractality can be observed from Figure 4.5.f (i.e., while the wideness of
the generalized Hurst exponent for the individual river flows is at least 1.5, the wide-
ness of generalized Hurst exponent from the multi-dimensional multi-fractal analysis
is narrower than 0.6). For completeness, we aim to investigate not only multiple river
flows exhibit similar dynamics, but also coupled river flows (measured across a com-
mon downstream path of water flow) that show different water levels (see Figure 4.5.a
and 4.5.b where the river flows is measured at several locations). From Figure 4.5.f
we notice that the generalized Hurst exponent for the measurements obtained across
the Sacramento river at Wilkins Slough CA while narrower than for the measurements
done at the other two sensor locations exhibits a more nonlinear shape. However, when
analyzing the scaling behavior of combined multiple time series the generalized Hurst
exponent shrinks and concentrates along a more symmetric step-like shape. In summary,
by employing the proposed mMFDFA method, we found that not only each of the river
flows have multi-fractal characteristics, but also that multiple coupled-locations across
river network basin can exhibit multi-fractal behavior.
4.5.2 Multi-fractality of water flow dynamics through a river basin
is season-dependent.
We investigate whether the multi-fractal behavior of the water flow through a large
river network changes over time (from season to season and from one year to another).
By applying the mMDFA method on multiple river flow time series obtained from the
Mississippi and Sacramento river networks and considering the water flow for several
years, we observe that the strength of multi-fractality (which is defined as the difference
between the maximum and minimum generalized Hurst exponents (discussed in Sec-
tion II)) changes over years and across seasons. For instance, Figure 3.a shows that the
72
Figure 4.5: a) Map of the California water resources including Sacramento River basins.
b) Map of Minnesota water resources including upper Mississippi River basins. c, d)
Time series of the water flow for three different measurement locations from the Missis-
sippi (c) and Sacramento (d) river basins in Oct 2007. e, f) Generalized Holder exponent
of individual river flows and combination of individual time series in part (c) and (d)
respectively.
73
Figure 4.6: a) The generalized Hurst exponent as a function of the q-th order moment for
multi-channel river stream data obtained from four different sites from the Mississippi
river basin during the month of September in different years, i.e., 2012, 2013 and 2014
respectively. b) The generalized Hurst exponent (encompassing the degree of multi-
fractality) as a function of the q-th order moment for multi-channel river stream data
obtained from four different sites within the Sacramento river in California during three
different months in 2014. c) The generalized Hurst exponent for the multi-channel river
stream data obtained from four different sites of Sacramento River in California during
the June-December time period in several consecutive years, i.e., 2012, 2013, and 2014.
shape and wideness of generalized Hurst exponent measured concomitantly for four sen-
sor locations across the Mississippi river basin changes significantly (e.g., the wideness
of the generalized Hurst exponent was X in 2012, Y in 2013 and Z in 2014). We show
the evolution of the multi-fractality of multi-variable river streams data from Mississippi
and Sacramento river basins in Figure 4.6.
4.5.3 Multi-fractal spectrum of the river network dynamics encom-
passes the effect of climate change.
We quantify the multi-fractal spectrums of several river flow time series. Although it
is intuitive that the river flow time series of a river stream at different locations and
branches should look similar, we observed that they are not the same and even the fluc-
tuation trends among them are different because of the effect of climate. We compare the
74
Figure 4.7: (a) And (b) Two sets of three time series representing the river flow levels
measured in the Mississippi river basin near Royalton, MN (Ch1), at St. Cloud, MN
(Ch2) and Mississippi River at HWY 610 in Brooklyn Park, MN (Ch3) (c) And (d). The
Generalized Hurst exponent as a function of the q-th order moment for the single- and
joint three-variable time series data obtained from the Mississippi River locations shown
in part (a) and (b) respectively.
multi-fractal spectrum of several river branches at different times. We show the multi-
fractal spectrums of multi-variable Mississippi River data at different times in Figure
4.7. Comparing Figure 4.7. a and Figure 4.7. b shows that as the fluctuation of differ-
ent river branches gets more similar, the multi-variable Generalized Hurst exponent gets
narrower and as the fluctuation of different locations gets less similar, the Generalized
Hurst exponent obtained by mMFDFA gets wider. The difference between fluctuations
of different sites of river stream is the result of environmental effects while having the
same trend comes from the fact that they belong to the same main river stream.
75
Chapter 5
Hardware prototyping of CPS
components
5.1 Hardware design challenges in MCPS
Advances in the physical sciences and engineering have enabled the development of
Cyber-Physical Systems (CPS). Health care systems require a special class of cyber-
physical systems that address their particular needs and constraints (Medical Cyber-
Physical Systems (MCPS)). Example of MCPSs includes artificial devices, monitoring
systems, robotic arms for surgeries, alarming systems, and other clinical or portable
devices. They are widely used both inside hospitals and also, they are used as portable
devices to improve the quality of life of patients. The effective use of MCPS reduces
the cost of health care by reducing the risk, improving the service delivery and reduc-
ing the human resources for delivering the health service. In previous Chapters, mostly
in Chapter 2, we have discussed the theory and design of the control algorithm for the
artificial pancreas. Though building the theory for control algorithm is the first step to
build a functional and reliable MCPS, hardware realization of them should also be con-
sidered closely To be able to produce those devices. We have classified the challenges
in hardware realization as follows:
Timing: Any feedback control algorithm should be able to perform the compu-
tation necessary for actuation in a timely manner. In other words, the control
algorithm should run fast enough to produce the necessary output in the feedback
76
control loop. To achieve the required timing , both sensing, control and actuation
should run on time. Otherwise, the device is not useful.
Energy: Energy and power consumption of medical devices is also a very critical
parameter. Many of these devices are battery operated and implantable. So, it
will not be practical if they will not be able to be operational for the time between
battery changing or charging periods. Also, for implantable devices, the tempera-
ture of the device should not exceed a limit because of the damage to the internal
organs and tissues. A lifetime of the devices is also affected by their power con-
sumption. In summary, the energy efficiency of medical devices is a critical design
parameter to take into account.
Device size: Because of being portable and implantable, sizing of the device is
also an important parameter. Medical devices should be able to fit the sizing
requirements for them to be practical.
In this Chapter, we bring the result of our prototype of the fractal model predictive
controller for the design of artificial pancreas.
5.2 Hardware prototyping of fractal optimal controllers
As described in Section 2.3, the optimal controller for blood glucose regulation can
be reduced to a linear program (LP). Bayliss et al. [105] investigated the hardware
implementation of an LP using the Simplex method. However, Simplex method does
not guarantee the convergence to the solution in polynomial time and so cannot be used
for real time control. In contrast, Interior-Point methods have fast polynomial-time
convergence which makes them suitable for real time control applications. Jerez et al
[15] investigated the hardware implementation of Interior-Point methods in the context
77
Figure 5.1: (a). Primal-dual formulation of interior point method. (b). Pseudo-code of
LU factorization. (c). Interior point algorithm. (d). Parallel architecture for updating
every matrix element.
of quadratic programming. Although applicable for real time control applications, this
implementation did not exploit the special format of matrices involved in the LP.
To avoid the Simplex method drawback and the Interior- Point inefficiency, we
developed an efficient hardware architecture of the proposed CPS-based AP capable
of exploiting the fractal characteristics and sparsity of the matrices involved in the LP.
For instance, by exploiting the characteristics of the fractional order derivatives and
appropriately ordering the variables within the interior point formulation, we reduced
the complexity of the matrix factorization step and eliminated the use of divide opera-
tions. These algorithmic approaches not only increase the speed of our implementation,
but also make our fractal optimal controller appropriate for hardware realization since
dividers can in- crease the hardware costs.
More precisely, our Interior-Point based hardware implementation of the fractal opti-
mal controller consists of a linear solver (i.e., solve Ax = b), a matrix to vector multi-
plication and an addition module, and a logic control module. The statement of primal
and dual problem for the interior point method is depicted in 5.1.a and the pseudo code
for the primal-dual interior point problem is depicted in 5.1.c. As can be seen there is
78
Figure 5.2: (a). Clock frequency and (b). number of resources on FPGA for the synthe-
sized controller and different problem sizes.
a linear system of equations which is repeatedly solved until the convergence criteria
are met. We have used the factorization algorithm depicted in 5.1.b and backward and
forward substitution for solving the linear system. Line 3 of the algorithm in 5.1.b,
requires division of all elements below diagonal by the corresponding diagonal element
in each column. However, by appropriate variable ordering there is no need for division
operation here and only negation of xi elements is enough and the rest of lower diagonal
elements remain unchanged. Elimination of divide operations in line 3 of 5.1.b not only
speedups the algorithm, but also makes this algorithm appropriate for hardware realiza-
tion since dividers are expensive in hardware. For line 7, which is the main part of the
factorization algorithm, we have used a parallel architecture which updates every matrix
element at each step by an adder and a multiplier depicted in 5.1.d. The proposed hard-
ware architecture was implemented in VHDL, prototyped on Xilinx Virtex-7 FPGAs
and the synthesized results are reported in 5.2.
79
5.3 Findings and suggestions for future work
In summary, by prototyping the fractional controller, we have shown that our proposed
theory works in practice as well. As future direction we propose investigating the how
the design parameters affects the performance of the controller:
Length of Control Horizon. We propose investigation of how the controller’s
features changes by changing the control Horizon. This reflects how it can adapt
to change in sensing and actuation technologies.
Change of Objective function. We propose investigating of alternative designs
for controller if we change the objective function and consider alternative objec-
tive functions.
Different technology. We have reported the result of hardware synthesis on
FPGA platform. We propose considering other platforms including micropro-
cessor.
80
Chapter 6
Conclusions and future research
directions
6.1 Major contributions of this Dissertation
Cyber-physical systems (CPS) rely on understanding the underlying physical compo-
nents to be used in engineered systems to finally built systems from synergy of com-
putational and physical components. Elucidating the dynamic nature of complex sys-
tems has important consequences not only on deciphering the topological interaction
among various system components, but also for predicting their function and behav-
ior. To date, most of the mathematical models for complex systems have ignored the
spatio-temporal (multi-scale) correlations and modeled their dynamics through memo-
ryless (integer order differential or difference) equations. Simply speaking, the current
mathematical models of complex systems neglect either the observed fractality in the
dynamics of individual components of the complex system or the higher-order cross-
correlations between various components of the complex system. To remedy this dis-
connect between the current mathematical modeling and the intrinsic characteristics of
complex systems, in this proposal, we address the modeling of the dynamics and then,
control of complex systems through fractal dimension analysis and fractional order con-
trol. We have also shown the advantages of our work through several applications espe-
cially from biological systems.
81
We consider the design of the control algorithm for artificial pancreas as a case
study to investigate the impact and importance of multifractal characteristics in
CPS design. We prove that the blood glucose (BG) dynamics is highly complex
and multifractal. We proposed a fractal-based model for BG in contrast to previ-
ous memory-less models and we show the advantage of our model in terms of the
goodness of fit. We formulate the fractional order model predictive controller in
artificial pancreas in contrast to conventional integer order controllers. We have
reported the advantage of use of this model in bringing the BG level to a reference
point even in the presence of uncertainty. We have also shown the feasibility of
implementation of fractal model predictive controller in hardware and we have
reported the hardware prototyping results.
We also extended our studies on gene expression networks as an example of a (bio-
logical) complex network. We investigate the scaling properties of gene expres-
sion time series in E. coli and S. cerevisiae. We investigate the individual gene
expression dynamics as well as the cross-dependency between them in the con-
text of gene regulatory network. Our results demonstrate that the gene expression
time series display fractal and long-range dependence characteristics. In addition,
the dynamics between genes and linked transcription factors in gene regulatory
networks are also fractal and long-range cross-correlated. The cross-correlation
exponents in gene regulatory networks are not unique. The distribution of the
cross-correlation exponents of gene regulatory networks for several types of cells
can be interpreted as a measure of the complexity of their functional behavior. We
extend the concept of cross-correlation fractal analysis and studied the dependen-
cies through concept of entropy of network by statistical distribution analysis on
cross-dependencies.
82
We also analyze the dynamical process in the complex networks. We investi-
gate the existence of fractional order dynamics in complex networks. We intro-
duce a multi- dimensional (-variable) multifractal detrended fluctuation analysis
(mMFDFa) method that allows quantifying the correlation complexity and mem-
ory existing between groups of interacting components of a complex system. Our
proposed method analyzes the higher-order statistics and quantifies the power
law scaling exhibited by dynamical fluctuations of coupled processes (time series
encompassing the evolution of complex system components). We have shown
how applying this method is useful on correlated and independent time series.
More precisely, our new method for computation of multi-fractal spectrum of
multi-variable time series can distinguish between correlate and independent time
series when put together for analysis. As an example of application this method,
we have reported our results for distinguishing different states of the human brain
in different sleep stages.
In summary, our contribution lies both on advancing theory and practice of the new
methods. We have formulated the model predictive fractal controller for application of
artificial pancreas and report the feasibility of implementation of this controller on real
hardware by prototyping. We have also, developed new modeling method for multi-
fractal detrended fluctuation of multi-variable time series. We have shown how this new
method is useful in deciphering complex interactions of components in the context of
complex networks. We have shown by extensive examples from both synthetic data and
dynamics of human brain through different sleep stages. We have also reported the frac-
tal nature of cross-dependencies in gene regulatory networks and show how extending
the notion of cross-correlation is useful in studying complex biological networks.
83
Figure 6.1: a. Singularity spectra of the of heart rate signals for healthy and dis-
eased hearts ([8]). b. Complexity index (sum of Sample Entropy values for scales 1
to 6, inclusively) and DFA-exponent (time scale greater than one hour) for the non-
diabetic and diabetic groups. Symbols with error bars represent mean and SD values,
respectively([9]).
6.2 Limitations and future research directions
In first part of this thesis, we have worked on formulating control algorithm for control
of blood glucose dynamics considering the fractal dynamics. We still have considered
the target value of the control signal (here blood glucose level) as the objective of the
optimization algorithm:
min
u(t)
t
f
R
0
(g(t)g
ref
(t))
2
dt (6.1)
Studies on several biological signals including heart rate and blood glucose dynam-
ics shows that healthy individuals have wider multi-fractal spectrum while the spectrum
for non-healthy individuals is narrow and is more like a mono-fractal signal. This moti-
vate us to rethink about control problem and think about the complexity of the biological
signals as the proffered output rather than only relying on target value of them. Hence,
complexity-aware control can be a new interesting direction for control algorithms:
84
6.2.1 Complexity Aware Control
Motivation. Studying physiological signals of healthy individuals reveals the more
complex nature of healthy vs. unhealthy statuses. This complexity is reported by multi-
fractal analysis and multi-scale entropy analysis [76] [62]. For example healthy indi-
viduals have multi- fractal heart rate dynamics while people suffering from heart dis-
eases have mono-fractal heart rate dynamics [62]. Another example is about dynamics
of blood glucose levels. Blood glucose dynamics of healthy individuals is reported to
be more complex than that of diabetes patience by multi-scale entropy analysis method
[77]. These facts have motivated us to look for new control objectives in medical devices
that do not only aim to bring the target physiological signals close to a reference value,
but also to maintain the signal complexity as well.
Problem Statement. Aiming to control a target state variable (for example blood
glucose level), we want to bring the state variable close to a reference value and also
maintain the complexity of the state variable. To do so, we need to bring new terms
in control objective in order to increase the complexity. In contrast to our formulation
for design of artificial pancreas shown in Figure 6.2 which only considers fractality in
dynamics of the blood glucose, the new approach should consider it in the objective
function as well. As a measure of the complexity, we can use the wideness of the
multi-fractal spectrum or number of Hurst exponents needed to describe the multi-fractal
spectrum. In order to formulate the complexity measure, we need to consider the time-
dependent fractional order derivative. Assuming that the state variable is brought close
to a reference value, we need to use another controller to keep the complexity of the
state variable. Given a finite control horizon, we want to maximize the number of time-
dependent fractional order derivative over the control horizon. We aim to formulate this
problem and investigate the convexity of the problem.
85
Figure 6.2: Our proposed model predictive fractal control formulation. We have
incorporate the fractal dynamics into constraints of the problem in constrast to classical
non-fractal optimal control problem.
Challenges. Investigating the computation complexity and convexity of the pro-
posed problem statement is challenging. We need to investigate the computation com-
plexity of the online estimation of the time-dependent fractional order deviate and fre-
quency of its repetition over the finite control horizon. This might affect how we want
to incorporate the complexity parameter in the objective function of the control algo-
rithm. Another important observation is how the new formulation affects the sampling
rate needed for feedback control. Since measuring multi-fractal spectrum requires to
have multiple data point over a time frame, the resolution of the multi-fractal spectrum
formulation in the objective function of the controller affects the sampling rate of mea-
surements and this can cause limitation on the performance of the method.
86
6.2.2 Controllability of Complex Networks
Motivation So far, we have addressed the control problem for a single state variable.
However, studying complex systems, which are a collection of a large number of compo-
nents, which interact with each other and the environment, requires considering complex
net- works. Examples of complex networks include biological, technological, ecolog-
ical or transportation networks. In complex networks, many components interact with
each other and the environment to self-organize into specific emergence structure. Mod-
eling dynamics of complex networks is the interest of several researches from several
communities. However, most of them have ignored the observed fractality in dynam-
ics of individual components and the collective behavior of network. Also, optimally
controlling the complex network in order to achieve certain tasks requires accurate math-
ematical models to include the observed fractality and long-range memory in complex
networks (e.g. biological networks).
Problem Statement Given N number of state variables corresponding to different
component of a complex network, the aim of the control problem is to bring the states
variables to target reference levels as close as possible in finite time horizon. In this
problem, different fractional order derivative corresponding to each state variable is
considered in the constraints of the optimal control problem. In summary, given con-
nectivity matrix A, the optimal control problem is the following:
Challenges One of the challenges in solving large-scale fractal optimization prob-
lem is computation complexity and scalability. We aim to investigate possible salable
hardware architectures to be able to reduce the run time of solving the optimal control
problem.
87
Reference List
[1] R. Schumer, M. M. Meerschaert, and B. Baeumer, “Fractional advection-
dispersion equations for modeling transport at the earth surface,”JournalofGeo-
physicalResearch: EarthSurface, vol. 114, no. F4, 2009.
[2] “Diabetes research in children network (direcnet). available online:
http://direcnet.jaeb.org.”
[3] M. Ghorbani and P. Bogdan, “A cyber-physical system approach to artificial pan-
creas design,” in2013InternationalConferenceonHardware/SoftwareCodesign
andSystemSynthesis(CODES+ISSS), pp. 1–10, IEEE, 2013.
[4] D. J. Lockhart and E. A. Winzeler, “Genomics, gene expression and dna arrays,”
Nature, vol. 405, no. 6788, p. 827, 2000.
[5] D. Marbach, J. C. Costello, R. K¨ uffner, N. M. Vega, R. J. Prill, D. M. Camacho,
K. R. Allison, A. Aderhold, R. Bonneau, Y . Chen, et al., “Wisdom of crowds for
robust gene network inference,”Naturemethods, vol. 9, no. 8, p. 796, 2012.
[6] H. Xiong and P. Shang, “Detrended fluctuation analysis of multivariate time
series,” Communications in Nonlinear Science and Numerical Simulation,
vol. 42, pp. 12–21, 2017.
[7] B. J. West, M. Latka, M. Glaubic-Latka, and D. Latka, “Multifractality of cere-
bral blood flow,”PhysicaA:StatisticalMechanicsanditsApplications, vol. 318,
no. 3-4, pp. 453–460, 2003.
[8] https://www.physionet.org/tutorials/multifractal/humanheart.htm.
[9] M. D. Costa, T. Henriques, M. N. Munshi, A. R. Segal, and A. L. Goldberger,
“Dynamical glucometry: use of multiscale entropy analysis in diabetes,” Chaos:
An Interdisciplinary Journal of Nonlinear Science, vol. 24, no. 3, p. 033139,
2014.
[10] M. M. Meerschaert and H.-P. Scheffler, “Semistable l´ evy motion,” Fractional
CalculusandAppliedAnalysis, vol. 5, no. 1, pp. 27–54, 2002.
88
[11] J. Radziuk, “The artificial pancreas,” Diabetes, vol. 61, no. 9, pp. 2221–2224,
2012.
[12] P. Eric Renard MD, J. Place, M. Cantwell, and C. C. Palerm, “Closed-loop insulin
delivery using a subcutaneous glucose sensor and intraperitoneal insulin deliv-
ery,”DiabetesCare, vol. 33, no. 1, p. 121, 2010.
[13] G. Marchetti, M. Barolo, L. Jovanovic, H. Zisser, and D. E. Seborg, “An improved
pid switching control strategy for type 1 diabetes,”ieeetransactionsonbiomedi-
calengineering, vol. 55, no. 3, pp. 857–865, 2008.
[14] L. Ljung, H. Hjalmarsson, and H. Ohlsson, “Four encounters with system identi-
fication,”EuropeanJournalofControl, vol. 17, no. 5-6, pp. 449–471, 2011.
[15] J. L. Jerez, G. A. Constantinides, and E. C. Kerrigan, “An fpga implementa-
tion of a sparse quadratic programming solver for constrained predictive control,”
in Proceedings of the 19th ACM/SIGDA international symposium on Field pro-
grammablegatearrays, pp. 209–218, ACM, 2011.
[16] B. B. Mandelbrot and M. M. Novak, Thinking in patterns: Fractals and related
phenomenainnature. World Scientific, 2004.
[17] G. M. Steil and J. Reifman, “Mathematical modeling research to support the
development of automated insulin-delivery systems,” 2009.
[18] P. Reichard, B.-Y . Nilsson, and U. Rosenqvist, “The effect of long-term intensi-
fied insulin treatment on the development of microvascular complications of dia-
betes mellitus,” New England Journal of Medicine, vol. 329, no. 5, pp. 304–309,
1993.
[19] M. O’Grady, P. John, and A. Winn, “Changes in medicare spending for type 1
diabetes with the introduction of the artificial pancreas,” New York (NY): JDRF,
2011.
[20] G. Sparacino, A. Facchinetti, and C. Cobelli, “smart continuous glucose monitor-
ing sensors: on-line signal processing issues,” Sensors, vol. 10, no. 7, pp. 6751–
6772, 2010.
[21] D. B. Keenan, B. Grosman, H. W. Clark, A. Roy, S. A. Weinzimer, R. V . Shah,
and J. J. Mastrototaro, “Continuous glucose monitoring considerations for the
development of a closed-loop artificial pancreas system,” 2011.
[22] C. Chen, Q. Xie, D. Yang, H. Xiao, Y . Fu, Y . Tan, and S. Yao, “Recent advances
in electrochemical glucose biosensors: a review,” Rsc Advances, vol. 3, no. 14,
pp. 4473–4491, 2013.
89
[23] L. Hoeks, W. Greven, and H. De Valk, “Real-time continuous glucose monitor-
ing system for treatment of diabetes: a systematic review,” Diabetic Medicine,
vol. 28, no. 4, pp. 386–394, 2011.
[24] J. Kildegaard, T. F. Christensen, and O. K. Hejlesen, “Sources of glycemic vari-
abilitywhat type of technology is needed?,”Journalofdiabetesscienceandtech-
nology, vol. 3, no. 4, pp. 986–991, 2009.
[25] A. Ouattara, A. Grimaldi, and B. Riou, “Blood glucose variability a new paradigm
in critical care?,” Anesthesiology: The Journal of the American Society of Anes-
thesiologists, vol. 105, no. 2, pp. 233–234, 2006.
[26] C. Cobelli, C. Dalla Man, G. Sparacino, L. Magni, G. De Nicolao, and B. P.
Kovatchev, “Diabetes: models, signals, and control,”IEEEreviewsinbiomedical
engineering, vol. 2, pp. 54–96, 2009.
[27] C. Ellingsen, E. Dassau, H. Zisser, B. Grosman, M. W. Percival, L. Jovanoviˇ c,
and F. J. Doyle III, “Safety constraints in an artificial pancreatic cell: an imple-
mentation of model predictive control with insulin on board,”Journalofdiabetes
scienceandtechnology, vol. 3, no. 3, pp. 536–544, 2009.
[28] F. H. El-Khatib, S. J. Russell, D. M. Nathan, R. G. Sutherlin, and E. R. Dami-
ano, “A bihormonal closed-loop artificial pancreas for type 1 diabetes,” Science
translationalmedicine, vol. 2, no. 27, pp. 27ra27–27ra27, 2010.
[29] D. Bruttomesso, A. Farret, S. Costa, M. C. Marescotti, M. Vettore, A. Avogaro,
A. Tiengo, C. Dalla Man, J. Place, A. Facchinetti, et al., “Closed-loop artificial
pancreas using subcutaneous glucose sensing and insulin delivery and a model
predictive control algorithm: preliminary studies in padova and montpellier,”
2009.
[30] D. R. Group et al., “Epidemiology of severe hypoglycemia in the diabetes con-
trol and complications trial,” The American journal of medicine, vol. 90, no. 4,
pp. 450–459, 1991.
[31] P. G. Jacobs, J. El Youssef, J. R. Castle, J. M. Engle, D. L. Branigan, P. Johnson,
R. Massoud, A. Kamath, and W. K. Ward, “Development of a fully automated
closed loop artificial pancreas control system with dual pump delivery of insulin
and glucagon,” in2011AnnualInternationalConferenceoftheIEEEEngineering
inMedicineandBiologySociety, pp. 397–400, IEEE, 2011.
[32] G. M. Steil, K. Rebrin, R. Janowski, C. Darwin, and M. F. Saad, “Modeling
-cell insulin secretion-implications for closed-loop glucose homeostasis,” Dia-
betestechnology&therapeutics, vol. 5, no. 6, pp. 953–964, 2003.
90
[33] J.-H. Cho, S. Sung, and I.-B. Lee, “Cascade control strategy for external carbon
dosage in predenitrifying process,” Water science and technology, vol. 45, no. 4-
5, pp. 53–60, 2002.
[34] A. Ali and S. Majhi, “Pid controller tuning for integrating processes,” ISA trans-
actions, vol. 49, no. 1, pp. 70–78, 2010.
[35] S. A. Weinzimer, G. M. Steil, K. L. Swan, J. Dziura, N. Kurtz, and W. V . Tambor-
lane, “Fully automated closed-loop insulin delivery versus semiautomated hybrid
control in pediatric patients with type 1 diabetes using an artificial pancreas,”
Diabetescare, vol. 31, no. 5, pp. 934–939, 2008.
[36] E. Renard, J. Place, M. Cantwell, H. Chevassus, and C. C. Palerm, “Closed-loop
insulin delivery using a subcutaneous glucose sensor and intraperitoneal insulin
delivery: feasibility study testing a new model for the artificial pancreas,” Dia-
betescare, vol. 33, no. 1, pp. 121–127, 2010.
[37] L. Magni, D. M. Raimondo, L. Bossi, C. Dalla Man, G. De Nicolao,
B. Kovatchev, and C. Cobelli, “Model predictive control of type 1 diabetes: an in
silico trial,” 2007.
[38] W. L. Clarke, S. Anderson, M. Breton, S. Patek, L. Kashmer, and B. Kovatchev,
“Closed-loop artificial pancreas using subcutaneous glucose sensing and insulin
delivery and a model predictive control algorithm: the virginia experience,” 2009.
[39] K. van Heusden, E. Dassau, H. C. Zisser, D. E. Seborg, and F. J. Doyle III,
“Control-relevant models for glucose control using a priori patient characteris-
tics,” IEEE transactions on biomedical engineering, vol. 59, no. 7, pp. 1839–
1849, 2012.
[40] H. Lee, B. A. Buckingham, D. M. Wilson, and B. W. Bequette, “A closed-loop
artificial pancreas using model predictive control and a sliding meal size estima-
tor,” 2009.
[41] R. Hovorka, J. Powrie, G. Smith, P. Sonksen, E. Carson, and R. Jones, “Five-
compartment model of insulin kinetics and its use to investigate action of
chloroquine in niddm,” American Journal of Physiology-Endocrinology and
Metabolism, vol. 265, no. 1, pp. E162–E175, 1993.
[42] K. Zarkogianni, A. Vazeou, S. G. Mougiakakou, A. Prountzou, and K. S. Nikita,
“An insulin infusion advisory system based on autotuning nonlinear model-
predictive control,”IEEEtransactionsonbiomedicalengineering, vol. 58, no. 9,
pp. 2467–2477, 2011.
91
[43] A. D. Association et al., “Accuracy of the glucowatch g2 biographer and the
continuous glucose monitoring system during hypoglycemia: experience of the
diabetes research in children network,”Diabetescare, vol. 27, no. 3, pp. 722–726,
2004.
[44] B. Whitcher and M. J. Jensen, “Wavelet estimation of a local long memory param-
eter,”ExplorationGeophysics, vol. 31, no. 1/2, pp. 94–103, 2000.
[45] “A.r. myers,medicine, lippincott williams & wilkins, medical, 2005..”
[46] S. A. Teichmann and M. M. Babu, “Gene regulatory network growth by duplica-
tion,”Naturegenetics, vol. 36, no. 5, p. 492, 2004.
[47] S. Huang, G. Eichler, Y . Bar-Yam, and D. E. Ingber, “Cell fates as high-
dimensional attractor states of a complex gene regulatory network,” Physical
reviewletters, vol. 94, no. 12, p. 128701, 2005.
[48] K. D¨ uvel, J. L. Yecies, S. Menon, P. Raman, A. I. Lipovsky, A. L. Souza, E. Tri-
antafellow, Q. Ma, R. Gorski, S. Cleaver, et al., “Activation of a metabolic gene
regulatory network downstream of mtor complex 1,” Molecular cell, vol. 39,
no. 2, pp. 171–183, 2010.
[49] R. Niedenthal, L. Riles, M. Johnston, and J. Hegemann, “Green fluorescent pro-
tein as a marker for gene expression and subcellular localization in budding
yeast,”Yeast, vol. 12, no. 8, pp. 773–786, 1996.
[50] R. Lanza, J. Gearhart, B. Hogan, D. Melton, R. Pedersen, E. D. Thomas, J. A.
Thomson, and M. West,Essentialsofstemcellbiology. Elsevier, 2005.
[51] T. Malone, R. Laubacher, and C. Dellarocas, “Mapping the genome of collec-
tive intelligence,” MIT Centre for Collective Intelligence Working Paper, vol. 1,
pp. 341–358, 2009.
[52] G. Kolesov, Z. Wunderlich, O. N. Laikova, M. S. Gelfand, and L. A. Mirny, “How
gene order is influenced by the biophysics of transcription regulation,” Proceed-
ings of the National Academy of Sciences, vol. 104, no. 35, pp. 13948–13953,
2007.
[53] O. Pulkkinen and R. Metzler, “Distance matters: the impact of gene proximity
in bacterial gene regulation,”Physicalreviewletters, vol. 110, no. 19, p. 198101,
2013.
[54] M. Bauer and R. Metzler, “In vivo facilitated diffusion model,” PloS one, vol. 8,
no. 1, p. e53956, 2013.
92
[55] J. Elf, G.-W. Li, and X. S. Xie, “Probing transcription factor dynamics at the
single-molecule level in a living cell,” Science, vol. 316, no. 5828, pp. 1191–
1194, 2007.
[56] T. E. Kuhlman and E. C. Cox, “Gene location and dna density determine tran-
scription factor distributions in escherichia coli,” Molecular systems biology,
vol. 8, no. 1, p. 610, 2012.
[57] M. Tsuchiya, A. Giuliani, M. Hashimoto, J. Erenpreisa, and K. Yoshikawa, “Self-
organizing global gene expression regulated through criticality: mechanism of the
cell-fate change,”PloSone, vol. 11, no. 12, p. e0167912, 2016.
[58] M. Tsuchyia, S. T. Wong, Z. X. Yeo, A. Colosimo, M. C. Palumbo, L. Farina,
M. Crescenzi, A. Mazzola, R. Negri, M. M. Bianchi, et al., “Gene expression
waves: cell cycle independent collective dynamics in cultured cells,” The FEBS
journal, vol. 274, no. 11, pp. 2878–2886, 2007.
[59] J. W. Kantelhardt, S. A. Zschiegner, E. Koscielny-Bunde, S. Havlin, A. Bunde,
and H. E. Stanley, “Multifractal detrended fluctuation analysis of nonstationary
time series,” Physica A: Statistical Mechanics and its Applications, vol. 316,
no. 1-4, pp. 87–114, 2002.
[60] J. W. Kantelhardt, S. A. Zschiegner, E. Koscielny-Bunde, S. Havlin, A. Bunde,
and H. E. Stanley, “Multifractal detrended fluctuation analysis of nonstationary
time series,” Physica A: Statistical Mechanics and its Applications, vol. 316,
no. 1-4, pp. 87–114, 2002.
[61] B. Podobnik and H. E. Stanley, “Detrended cross-correlation analysis: a new
method for analyzing two nonstationary time series,” Physical review letters,
vol. 100, no. 8, p. 084102, 2008.
[62] B. B. Mandelbrot, A. J. Fisher, and L. E. Calvet, “A multifractal model of asset
returns,” 1997.
[63] A. Arneodo, E. Bacry, and J.-F. Muzy, “Random cascades on wavelet dyadic
trees,”JournalofMathematicalPhysics, vol. 39, no. 8, pp. 4142–4164, 1998.
[64] C. E. Shannon, “A mathematical theory of communication. acm sigmobile mob,”
Comput.Commun.Rev, vol. 5, no. 1, pp. 3–55, 2001.
[65] K. Anand and G. Bianconi, “Entropy measures for networks: Toward an informa-
tion theory of complex topologies,”PhysicalReviewE, vol. 80, no. 4, p. 045102,
2009.
93
[66] K. Anand, G. Bianconi, and S. Severini, “Shannon and von neumann entropy
of random networks with heterogeneous expected degree,” Physical Review E,
vol. 83, no. 3, p. 036109, 2011.
[67] J. W. Kantelhardt, E. Koscielny-Bunde, H. H. Rego, S. Havlin, and A. Bunde,
“Detecting long-range correlations with detrended fluctuation analysis,” Physica
A: Statistical Mechanics and its Applications, vol. 295, no. 3-4, pp. 441–454,
2001.
[68] O. Muzychuk, “Using bimodal probability distributions in the problems of brow-
nian diffusion,” Radiophysics and quantum electronics, vol. 49, no. 8, pp. 645–
655, 2006.
[69] B. Efron,Thejackknife,thebootstrap,andotherresamplingplans, vol. 38. Siam,
1982.
[70] C.-K. Peng, S. V . Buldyrev, S. Havlin, M. Simons, H. E. Stanley, and A. L.
Goldberger, “Mosaic organization of dna nucleotides,”Physicalreviewe, vol. 49,
no. 2, p. 1685, 1994.
[71] A. Arneodo, Y . d’Aubenton Carafa, E. Bacry, P. Graves, J. Muzy, and C. Ther-
mes, “Wavelet based fractal analysis of dna sequences,” Physica D: Nonlinear
Phenomena, vol. 96, no. 1-4, pp. 291–320, 1996.
[72] Q. Zhang, S. Zhou, and X. Wei, “An efficient approach for dna fractal-based
image encryption,” Applied Mathematics & Information Sciences, vol. 5, no. 3,
pp. 445–459, 2011.
[73] K. Kiyono, Z. R. Struzik, and Y . Yamamoto, “Estimator of a non-gaussian param-
eter in multiplicative log-normal models,” Physical Review E, vol. 76, no. 4,
p. 041113, 2007.
[74] E. Bacry and J. F. Muzy, “Log-infinitely divisible multifractal processes,” Com-
municationsinMathematicalPhysics, vol. 236, no. 3, pp. 449–475, 2003.
[75] P. Chainais, R. Riedi, and P. Abry, “On non-scale-invariant infinitely divisible
cascades,” IEEE transactions on Information Theory, vol. 51, no. 3, pp. 1063–
1083, 2005.
[76] Y . Xue and P. Bogdan, “Constructing compact causal mathematical models
for complex dynamics,” in Proceedings of the 8th International Conference on
Cyber-PhysicalSystems, pp. 97–107, ACM, 2017.
[77] A. E. Teschendorff and T. Enver, “Single-cell entropy for accurate estimation
of differentiation potency from a cells transcriptome,” Nature communications,
vol. 8, p. 15599, 2017.
94
[78] N. Ji, T. C. Middelkoop, R. A. Mentink, M. C. Betist, S. Tonegawa, D. Mooijman,
H. C. Korswagen, and A. van Oudenaarden, “Feedback control of gene expres-
sion variability in the caenorhabditis elegans wnt pathway,” Cell, vol. 155, no. 4,
pp. 869–880, 2013.
[79] W. Maass, “Noise as a resource for computation and learning in networks of
spiking neurons,”ProceedingsoftheIEEE, vol. 102, no. 5, pp. 860–880, 2014.
[80] P. Bernaola-Galv´ an, R. Rom´ an-Rold´ an, and J. L. Oliver, “Compositional segmen-
tation and long-range fractal correlations in dna sequences,” Physical Review E,
vol. 53, no. 5, p. 5181, 1996.
[81] G. B. West, J. H. Brown, and B. J. Enquist, “The origin of universal scaling laws
in biology,”Scalinginbiology, pp. 87–112, 2000.
[82] G. B. West, J. Brown, and B. Enquist, “Scaling in biology: patterns and processes,
causes and consequences,”Scalinginbiology, pp. 87–112, 2000.
[83] Y .-Y . Liu, J.-J. Slotine, and A.-L. Barab´ asi, “Controllability of complex net-
works,”nature, vol. 473, no. 7346, p. 167, 2011.
[84] R. Albert and A.-L. Barab´ asi, “Dynamics of complex systems: Scaling laws for
the period of boolean networks,”PhysicalReviewLetters, vol. 84, no. 24, p. 5660,
2000.
[85] M. A. de Menezes and A.-L. Barab´ asi, “Separating internal and external dynam-
ics of complex systems,”Physicalreviewletters, vol. 93, no. 6, p. 068701, 2004.
[86] R. Elber, A. Ghosh, and A. Cardenas, “Long time dynamics of complex systems,”
Accountsofchemicalresearch, vol. 35, no. 6, pp. 396–403, 2002.
[87] P. C. Ivanov, L. A. N. Amaral, A. L. Goldberger, S. Havlin, M. G. Rosenblum,
Z. R. Struzik, and H. E. Stanley, “Multifractality in human heartbeat dynamics,”
Nature, vol. 399, no. 6735, p. 461, 1999.
[88] B. Podobnik and H. E. Stanley, “Detrended cross-correlation analysis: a new
method for analyzing two nonstationary time series,” Physical review letters,
vol. 100, no. 8, p. 084102, 2008.
[89] M. Ghorbani and P. Bogdan, “A cyber-physical system approach to artificial pan-
creas design,” in2013InternationalConferenceonHardware/SoftwareCodesign
andSystemSynthesis(CODES+ISSS), pp. 1–10, IEEE, 2013.
[90] S. Makeig, S. Debener, J. Onton, and A. Delorme, “Mining event-related brain
dynamics,”Trendsincognitivesciences, vol. 8, no. 5, pp. 204–210, 2004.
95
[91] T.-P. Jung, S. Makeig, M. J. McKeown, A. J. Bell, T.-W. Lee, T. J. Sejnowski,
et al., “Imaging brain dynamics using independent component analysis,” Pro-
ceedingsoftheIEEE, vol. 89, no. 7, pp. 1107–1122, 2001.
[92] A.-M. Wink, E. Bullmore, A. Barnes, F. Bernard, and J. Suckling, “Monofrac-
tal and multifractal dynamics of low frequency endogenous brain oscillations in
functional mri,”Humanbrainmapping, vol. 29, no. 7, pp. 791–801, 2008.
[93] J. Suckling, A. M. Wink, F. A. Bernard, A. Barnes, and E. Bullmore, “Endoge-
nous multifractal brain dynamics are modulated by age, cholinergic blockade
and cognitive performance,” Journal of neuroscience methods, vol. 174, no. 2,
pp. 292–300, 2008.
[94] B. Podobnik and H. E. Stanley, “Detrended cross-correlation analysis: a new
method for analyzing two nonstationary time series,” Physical review letters,
vol. 100, no. 8, p. 084102, 2008.
[95] X. Zhang, M. Zeng, and Q. Meng, “Multivariate multifractal detrended fluctua-
tion analysis of 3d wind field signals,” Physica A: Statistical Mechanics and its
Applications, vol. 490, pp. 513–523, 2018.
[96] B. Podobnik, D. Horvatic, A. L. Ng, H. E. Stanley, and P. C. Ivanov, “Model-
ing long-range cross-correlations in two-component arfima and fiarch processes,”
PhysicaA:StatisticalMechanicsanditsApplications, vol. 387, no. 15, pp. 3954–
3959, 2008.
[97] L. K. Gallos, H. A. Makse, and M. Sigman, “A small world of weak ties provides
optimal global integration of self-similar modules in functional brain networks,”
ProceedingsoftheNationalAcademyofSciences, vol. 109, no. 8, pp. 2825–2830,
2012.
[98] L. Gallos, M. Sigman, and H. Makse, “The conundrum of functional brain net-
works: small-world efficiency or fractal modularity,” Frontiers in physiology,
vol. 3, p. 123, 2012.
[99] A. Clark,Beingthere: Puttingbrain,body,andworldtogetheragain. MIT press,
1998.
[100] R. Lande, “Quantitative genetic analysis of multivariate evolution, applied to
brain: body size allometry,”Evolution, vol. 33, no. 1Part2, pp. 402–416, 1979.
[101] K. J. Maloney, E. Cape, J. Gotman, and B. Jones, “High-frequency
elec-
troencephalogram activity in association with sleep-wake states and spontaneous
behaviors in the rat,”Neuroscience, vol. 76, no. 2, pp. 541–555, 1997.
96
[102] J.-E. Kang, M. M. Lim, R. J. Bateman, J. J. Lee, L. P. Smyth, J. R. Cirrito,
N. Fujiki, S. Nishino, and D. M. Holtzman, “Amyloid- dynamics are regulated
by orexin and the sleep-wake cycle,”Science, vol. 326, no. 5955, pp. 1005–1007,
2009.
[103] D. C. Piper, N. Upton, M. I. Smith, and A. J. Hunter, “The novel brain neu-
ropeptide, orexin-a, modulates the sleep–wake cycle of rats,” European Journal
ofNeuroscience, vol. 12, no. 2, pp. 726–730, 2000.
[104] M. G. Terzano, L. Parrino, A. Smerieri, R. Chervin, S. Chokroverty, C. Guillem-
inault, M. Hirshkowitz, M. Mahowald, H. Moldofsky, A. Rosa, et al., “Atlas,
rules, and recording techniques for the scoring of cyclic alternating pattern (cap)
in human sleep,”Sleepmedicine, vol. 3, no. 2, pp. 187–199, 2002.
[105] S. Bayliss, G. A. Constantinides, W. Luk, et al., “An fpga implementation of
the simplex algorithm,” in 2006 IEEE International Conference on Field Pro-
grammableTechnology, pp. 49–56, IEEE, 2006.
97
Abstract (if available)
Abstract
Elucidating the dynamic nature of complex systems has significant consequences not only on deciphering the topological interaction among various system components but also for predicting their function and behavior. To date, most of the mathematical models for complex systems have ignored the spatio-temporal (multi-scale) correlations and modeled their dynamics through memoryless (integer order differential or difference) equations. Directly speaking, the current mathematical models of complex systems neglect either the observed fractality in the dynamics of individual components of the complex system or the higher-order cross-correlations between various components of the complex system. To remedy this disconnect between the current mathematical modeling and the intrinsic characteristics of complex systems, in this thesis, we address the modeling of the dynamics and then, control of complex systems through fractal dimension analysis and fractional order control. To investigate the impact and importance of multi-fractal characteristics in CPS design, we consider the design of the control algorithm for artificial pancreas as a case study. We prove that the blood glucose (BG) dynamics are highly complex and multi-fractal. We proposed a fractal-based model for BG in contrast to previous memory-less models, and we show the advantage of our model in terms of the goodness of fit. We formulate the fractional order model predictive controller in the artificial pancreas in contrast to conventional integer order controllers. ❧ We have reported the advantage of the use of this model in bringing the BG level to a reference point even in the presence of uncertainty. We have also shown the feasibility of implementation of the fractal model predictive controller in hardware, and we have reported the hardware prototyping results. We also extended our studies on gene expression networks as an example of a (biological) complex network. We investigate the scaling properties of gene expression time series in E. coli and S. cerevisiae. We investigate the individual gene expression dynamics as well as the cross-dependency between them in the context of gene regulatory network. Our results demonstrate that the gene expression time series display fractal and long-range dependence characteristics. Besides, the dynamics between genes and linked transcription factors in gene regulatory networks are also fractal and long-range cross-correlated. The cross-correlation exponents in gene regulatory networks are not unique. The distribution of the cross-correlation exponents of gene regulatory networks for several types of cells can be interpreted as a measure of the complexity of their functional behavior. We extend the concept of cross-correlation fractal analysis and studied the dependencies through the concept of network entropy by statistical distribution analysis on cross-dependencies. We also analyze the dynamical process in complex networks. We investigate the existence of fractional order dynamics in complex networks. We introduce a multidimensional (-variable) multi-fractal detrended fluctuation analysis (mMFDFa) method that allows quantifying the correlation complexity and memory exhibited between groups of interacting components of a complex system. Our proposed method analyzes the higher-order statistics and quantifies the power law scaling exhibited by dynamical fluctuations of coupled processes (time series encompassing the evolution of complex system components). To illustrate the depth and benefits of this mathematical method, we investigate the multi-fractal spatio-temporal characteristics. First, we show how using this method is beneficial for synthetic correlated and uncorrelated fractal signals. Then, we used this analysis to distinguish the activity of human brain consisting of groups of neurons during various sleep periods. ❧ As future direction toward this thesis, we propose to address a few new problems. First, we want to formulate the complexity-aware control algorithm. Motivated by the observed complex multi-fractal dynamics in physiological signals, there can be a new class of control algorithms that can maintain the complexity of the target control variable(s). Second, addressing the problem of controllability of complex networks considering fractal dynamics will also have a significant contribution to the field. Lastly, investigating the scalability of the proposed fractal controllers and their hardware implementation is also crucial to make the theory into practice.
Linked assets
University of Southern California Dissertations and Theses
Conceptually similar
PDF
Theoretical and computational foundations for cyber‐physical systems design
PDF
Theoretical foundations and design methodologies for cyber-neural systems
PDF
Theoretical foundations for modeling, analysis and optimization of cyber-physical-human systems
PDF
Verification, learning and control in cyber-physical systems
PDF
Deterministic evolutionary game theory models for the nonlinear dynamics and control of chemotherapeutic resistance
PDF
Biological geometry-aware deep learning frameworks for enhancing medical cyber-physical systems
PDF
Analysis, design, and optimization of large-scale networks of dynamical systems
PDF
Data-driven aggregated fractal dynamic load modeling for early warning of voltage collapse
PDF
Dynamic graph analytics for cyber systems security applications
PDF
Towards green communications: energy efficient solutions for the next generation cellular mobile communication systems
PDF
Dealing with unknown unknowns
PDF
Defending industrial control systems: an end-to-end approach for managing cyber-physical risk
PDF
Data-driven and logic-based analysis of learning-enabled cyber-physical systems
PDF
Molecular dynamics studies of protein aggregation in unbounded and confined media
PDF
Stochastic dynamic power and thermal management techniques for multicore systems
PDF
Learning and decision making in networked systems
PDF
Architectures and algorithms of charge management and thermal control for energy storage systems and mobile devices
PDF
Acceleration of deep reinforcement learning: efficient algorithms and hardware mapping
PDF
Stochastic and multiscale models for urban and natural ecology
PDF
High-accuracy adaptive vibrational control for uncertain systems
Asset Metadata
Creator
Ghorbani, Mahboobeh
(author)
Core Title
Understanding dynamics of cyber-physical systems: mathematical models, control algorithms and hardware incarnations
School
Viterbi School of Engineering
Degree
Doctor of Philosophy
Degree Program
Electrical Engineering
Publication Date
08/15/2019
Defense Date
08/15/2019
Publisher
University of Southern California
(original),
University of Southern California. Libraries
(digital)
Tag
control algorithms,cyber-physical systems,dynamics,hardware incarnations,Models,OAI-PMH Harvest
Format
application/pdf
(imt)
Language
English
Contributor
Electronically uploaded by the author
(provenance)
Advisor
Bogdan, Paul (
committee chair
), Gupta, Sandeep (
committee member
), Jonckheere, Edmond (
committee member
), Nakano, Aiichiro (
committee member
)
Creator Email
mahboobe.ghorbani@gmail.com,mahboobg@usc.edu
Permanent Link (DOI)
https://doi.org/10.25549/usctheses-c89-211680
Unique identifier
UC11662796
Identifier
etd-GhorbaniMa-7800.pdf (filename),usctheses-c89-211680 (legacy record id)
Legacy Identifier
etd-GhorbaniMa-7800.pdf
Dmrecord
211680
Document Type
Dissertation
Format
application/pdf (imt)
Rights
Ghorbani, Mahboobeh
Type
texts
Source
University of Southern California
(contributing entity),
University of Southern California Dissertations and Theses
(collection)
Access Conditions
The author retains rights to his/her dissertation, thesis or other graduate work according to U.S. copyright law. Electronic access is being provided by the USC Libraries in agreement with the a...
Repository Name
University of Southern California Digital Library
Repository Location
USC Digital Library, University of Southern California, University Park Campus MC 2810, 3434 South Grand Avenue, 2nd Floor, Los Angeles, California 90089-2810, USA
Tags
control algorithms
cyber-physical systems
dynamics
hardware incarnations