# Exploring the Interplay of Intrinsic Fluctuation and Complexity in
Intracellular Calcium Dynamics

Athokpam Langlen Chanu<EMAIL_ADDRESS>(Corresponding author) Asia
Pacific Center for Theoretical Physics, Pohang, 37673, Republic of Korea  R.
K. Brojen Singh<EMAIL_ADDRESS>School of Computational and Integrative
Sciences, Jawaharlal Nehru University, New Delhi, 110067, India  Jae-Hyung
Jeon<EMAIL_ADDRESS>Asia Pacific Center for Theoretical Physics,
Pohang, 37673, Republic of Korea  Department of Physics, Pohang University of
Science and Technology (POSTECH), Pohang, 37673, Republic of Korea

###### Abstract

The concentration of intracellular calcium ion (Ca2+) exhibits complex
oscillations, including bursting and chaos, as observed experimentally. These
dynamics are influenced by inherent fluctuations within cells, which serve as
crucial determinants in cellular decision-making processes and fate
determination. In this study, we systematically explore the interplay between
intrinsic fluctuation and the complexity of intracellular cytosolic Ca2+
dynamics. To investigate this interplay, we employ complexity measures such as
permutation entropy and statistical complexity. Using the chemical Langevin
equation, we simulate the stochastic dynamics of cytosolic Ca2+. Our findings
reveal that permutation entropy and statistical complexity effectively
characterize the diverse, dynamic states of cytosolic Ca2+ and illustrate
their interactions with intrinsic fluctuation. Permutation entropy analysis
elucidates that the chaotic state is more sensitive to intrinsic fluctuation
than the other periodic states. Furthermore, we identify distinct states of
cytosolic Ca2+ occupying specific locations within the theoretical bounds of
the complexity-entropy causality plane. These locations indicate varying
complexity and information content as intrinsic fluctuation varies. When
adjusting the permutation order, the statistical complexity measure for the
different periodic and chaotic states exhibits peaks in an intermediate range
of intrinsic fluctuation values. Additionally, we identify scale-free or self-
similar patterns in this intermediate range, which are further corroborated by
multifractal detrended fluctuation analysis. These high-complexity states
likely correspond to optimal Ca2+ dynamics with biological significance,
revealing rich and complex dynamics shaped by the interplay of intrinsic
fluctuation and complexity. Our investigation enhances our understanding of
how intrinsic fluctuation modulates the complexity of intracellular Ca2+
dynamics that play crucial roles in biological cells.

Complex systems, complexity, calcium oscillations, chemical Langevin equation,
permutation entropy, complexity-entropy causality plane

##  I Introduction

Complex systems comprise numerous interacting sub-units or components,
characterized by intricate mutual interactions [1]. These interactions give
rise to non-linear phenomena, including bifurcation [2], which is considered
an elementary act of complexity [3, 4]. Other intriguing non-linear phenomena
encompass multistability [5], chaos [6], and fractals [7]. At microscopic and
mesoscopic scales, diverse biological systems showcase fluctuation-driven
stochastic dynamics, observed in biological cells [8], neural firing [9],
biochemical networks [10], genetic oscillator networks [11], molecular motors
[12], and more. Biological cells, complex biochemical systems [13, 14],
exchange energy and matter with the surrounding environment, operating far
from thermodynamic equilibrium. They exhibit dissipative structures [15], such
as temporal oscillations, waves, and patterns [16, 17, 18], which exhibit a
blend of orderliness (coherence) and randomness [19]. Exploring the complexity
of dissipative structures like biochemical oscillations is intriguing and
fundamentally significant.

In complex biochemical systems like living cells, fluctuations manifest as
intrinsic and extrinsic fluctuations [20]. Extrinsic fluctuations result from
external factors, whereas intrinsic fluctuations arise from inherent random
molecular interactions within a chemically reacting system [21]. These
interactions correspond to finite biochemical reactions within each
biochemical pathway of the cellular circuit, forming the basis for
intermolecular cross-talks that govern various cellular functions [22]. The
intrinsic fluctuations stemming from random molecular interactions play a
pivotal role in regulating cellular organizations [23]. These intrinsic
fluctuations can be modulated by the system size [24] (denoted as $V$), with
intrinsic fluctuation scaling approximately as $\sim\frac{1}{\sqrt{V}}$ [25].
Biological cells can undergo significant changes in size, ranging from a
$10$–$50\%$ expansion or shrinkage from their original size. Examples include
hippocampal neurons $\sim 10$–$45\%$ [26], erythrocytes (bone marrow cells)
$\sim 14$–$20\%$ [27], and somatic cells $\sim 20$–$40\%$ [28]. Intrinsic
fluctuations play a crucial role in cellular decision-making and fate
determination [29, 30]. Also, the collective behavior of dynamic structures,
encompassing stable, steady states, oscillations, chaos, etc., constitutes the
foundation of emergent molecular phenotype [31], shaping the overall
phenotypic nature of an organism [32]. Apart from simple periodic
oscillations, biological systems often exhibit complex oscillations like
bursting, birhythmicity, multi-periodicity, quasi-periodicity, and chaos.
These phenomena may correspond to complex molecular cross-talks and
distributions [33]. However, the relation between such complex dynamics and
molecular cross-talks remains not yet fully understood. In this regard, the
impact of intrinsic fluctuations on the complexity of dynamic states poses a
fundamental question crucial for gaining deeper insights into the complex
mechanisms orchestrated by biological systems. We aim to address this
fundamental question through a numerical study of complex oscillations and
chaos in a non-linear model of intracellular calcium ion (Ca2+) oscillation
based on the Ca2+-induced Ca2+ release (CICR) mechanism proposed by Houart et
al [34].

Calcium ion (Ca2+) functions as a vital messenger within biological cells,
undergoing intracellular Ca2+ oscillations [35] in various cell types,
including pancreatic cells [36, 37], hepatocytes [38], muscle cells [39, 40],
and neurons [41]. These oscillations play crucial roles not only in signal
transduction inside the cell [42, 43] but also in various physiological
processes, including gene expression [44], cell proliferation [45], and
neuronal differentiation [46]. The system of intracellular Ca2+ oscillations
is a non-equilibrium system [47]. Mathematical models have been developed to
explain the experimentally observed complex patterns of intracellular Ca2+
oscillations. In particular, Houart et al. [34] has developed a model based on
the non-linear feedback process of the Ca2+-induced Ca2+ release (CICR)
mechanism [48], prevalent in various cell types, including hepatocytes [49]
and cardiac [50] cells. In the CICR mechanism, the release of Ca2+ from
intracellular stores into the cytosol is activated by inositol trisphosphate
(InsP3) and cytosolic Ca2+ itself. This autocatalytic process produces a
variety of complex dynamical behaviors in the temporal patterns of Ca2+
oscillation. Intrinsic stochasticity in intracellular Ca2+ oscillations arises
due to finite cell size and a small number of reactants [51]. The interplay of
external noise strength and oscillatory dynamics has been studied in the
adaption of the biological system of Physarum polycephalum [52]. However, to
our knowledge, a systematic analysis of the interplay between intrinsic
fluctuation and the “complexity” of dynamical behavior such as complex
oscillations, including bursting, multi-periodicity, quasi-periodicity, and
chaos in the intracellular Ca2+ oscillations, has not been addressed earlier.
Complexity is a multifaceted notion [53, 54], with various complexity measures
proposed across different disciplines, including algorithmic complexity,
effective complexity, statistical complexity, fractal dimension, entropy,
degree of organization [55], among others. For a deterministic dynamical
system described by coupled, non-linear ordinary differential equations and
undergoing the bifurcation phenomenon, complexity is generally characterized
using Lyapunov exponents [2]. In information theory, the degree of disorder is
quantified using entropy [56]. Connecting non-linear dynamics and information
theory, Bandt and Pompe have proposed permutation entropy $H$ based on Shannon
entropy [57] as a measure of complexity in non-linear time series [58].
L$\mathrm{\acute{o}}$pez-Ruiz et al. [59] have defined a statistical measure
of complexity, denoted by $C$, relating order and information. A diagram of
$C$ versus $H$, where $H$ is regarded as an arrow of time, is known as the
complexity-entropy (CH) causality plane [60]. In this study, we use
permutation entropy $H$ and statistical complexity $C$ as quantitative
measures of complexity to analyze the complexities of various dynamical
behaviors of intracellular Ca2+ and their interplay with intrinsic
fluctuations.

We simulate the dynamics of the intracellular Ca2+ oscillation model proposed
by Houart et al. [34], employing a stochastic approach based on the chemical
Langevin equation (CLE) [25] to investigate the behavior of the stochastic
Ca2+ dynamics driven by intrinsic fluctuations. To delve into the intricate
relationship between intrinsic fluctuations and the complexity of
intracellular Ca2+ dynamics, we use permutation entropy [58] and the
statistical complexity measure [59]. We show that intrinsic fluctuations play
a significant role in modulating the patterns of Ca2+ oscillations within a
cell, even to the extent of disrupting the coherence of these oscillations.
Applying permutation entropy finds that at large intrinsic fluctuations,
states become noisier, and the distinction between the periodic and chaotic
states diminishes. Permutation entropy analysis elucidates that the chaotic
state is more sensitive to intrinsic fluctuation than the other periodic
states. Additionally, we identify distinct periodic and chaotic states of the
cytosolic Ca2+ at different planar locations on the CH causality plane as
intrinsic fluctuation varies, indicating varying complexity and information
content. We also observe peaks in the statistical complexities of the
different states within an intermediate range of values of intrinsic
fluctuation (or system size). Within this regime where heterogeneous
structures appear, we uncover scale-free or self-similar patterns in the
temporal oscillatory structures of cytosolic Ca2+. Our results demonstrate the
intricate interplay of intrinsic fluctuation and complexity in the periodic
and chaotic states of intracellular Ca2+ dynamics, signifying the regulatory
role of intrinsic fluctuations within a cell.

The structure of the paper unfolds as follows. In Sec. II, we describe the
model of intracellular Ca2+ oscillations adopted for this study. We outline
the methodologies employed, namely, the chemical Langevin equation (CLE) in
Sec. III.1, the permutation entropy (PE) in Sec. III.2, and the complexity-
entropy (CH) causality plane in Sec. III.3. Sec. IV presents and discusses our
main results. We then summarize our results with conclusions in Sec. V.

Figure 1: Schematic representation of intracellular calcium (Ca2+)
oscillations based on the mechanism of Ca2+-induced Ca2+-release (CICR)
(figure adapted from Ref. [49]).

##  II Model: Intracellular Calcium Oscillations

Consider a cell of system size $V$ (see Fig. 1). The variables $X,Y$, and $Z$
represent the population of free Ca2+ in the cytosol (cytosolic Ca2+), Ca2+
stored in the internal pool, and inositol 1,4,5-trisphosphate (InsP3),
respectively. Suppose $x=X/V,\ y=Y/V,$ and $z=Z/V$ represent the
concentrations of cytosolic Ca2+, stored Ca2+ and InsP3, respectively. The
time-evolution of these concentrations is governed by the following set of
coupled, non-linear ordinary differential equations (ODEs) [34]:

$\displaystyle\begin{aligned} \frac{dx}{dt}&=V_{0}+V_{1}\beta-
V_{2}+V_{3}+k_{f}y-kx,\\\ \frac{dy}{dt}&=V_{2}-V_{3}-k_{f}y,\\\
\frac{dz}{dt}&=\beta V_{4}-V_{5}-\epsilon z,\end{aligned}$ (1)

where

$\displaystyle V_{2}$ $\displaystyle=V_{M2}\frac{x^{2}}{k_{2}^{2}+x^{2}},$
$\displaystyle V_{3}$
$\displaystyle=V_{M3}\frac{x^{m}}{k_{x}^{m}+x^{m}}\frac{y^{2}}{k_{y}^{2}+y^{2}}\frac{z^{4}}{k_{z}^{4}+z^{4}},$
$\displaystyle V_{5}$
$\displaystyle=V_{M5}\frac{z^{p}}{k_{5}^{p}+z^{p}}\frac{x^{n}}{k_{d}^{n}+x^{n}}.$

In the above set of ODEs (1), $V_{0}$ represents the constant Ca2+ supply from
the extracellular medium, and the parameter $\beta$ denotes the degree of cell
stimulation by an agonist (e.g., hormone or neurotransmitter). The rate
$V_{2}\ (V_{3})$ represents the pumping of Ca2+ from the cytosol into the
internal pool (release of Ca2+ from the internal pool to the cytosol).
$V_{M2}$ and $V_{M3}$ are the maximum values of these rates. The parameters
$k_{2},\ k_{y},\ k_{x}$, and $k_{z}$ denote the threshold values for pumping,
release, and activation of release by Ca2+ and InsP3, respectively. The rate
constant $k_{f}$ measures the passive, linear leak of $y$ into $x$, and $k$
denotes the linear transport of cytosolic Ca2+ into the extracellular medium.
$V_{4}$ represents the maximum rate of stimulus-induced synthesis of InsP3,
$V_{5}$ represents the phosphorylation rate of InsP3 by the 3-kinase. $k_{5}$
denotes the half-saturation constant, and $k_{d}$ corresponds to a threshold
Ca2+ level. The parameters $m,n$, and $p$ are the Hill coefficients related to
cooperative processes. By adjusting the values of the rate constants and other
parameters involved in the model (1), the concentrations $x$, $y$, and $z$
exhibit various kinds of dynamics, including stable steady-state, simple
periodic oscillation, complex oscillations such as bursting, period doubling,
quasi-periodicity, and chaos [34].

##  III Methods

###  III.1 Stochastic Modeling with Chemical Langevin Equation (CLE)

Consider a well-stirred chemically reacting system of size $V$ comprising
three molecular species $X$, $Y$, and $Z$. Suppose the system is maintained at
a constant temperature $T$. The state vector of the populations of $X$, $Y$,
and $Z$ is denoted by
$\mathcal{S}=\mathcal{S}(t)=[X(t),Y(t),Z(t)]^{\mathrm{T}}$, where T denotes
transpose. The molecular interactions of $X$, $Y$, and $Z$ are modeled through
a set of $M$ chemical reaction channels:
$\ce{$a_{1}X$+$b_{1}Y$+$c_{1}Z$->[$k_{j}$]$a_{2}X$+$b_{2}Y$+$c_{2}Z$}$, where
$k_{j}\ ;\ j=1,2,\dots,M$ denotes the classical rate constant of the
$j^{\mathrm{th}}$ reaction. The sets $\\{a_{1},b_{1},c_{1}\\}$ and
$\\{a_{2},b_{2},c_{2}\\}$ represent the number of reactant and product
molecules, respectively. The propensity function $a_{j}$ is defined as
$a_{j}=c_{j}h_{j}$ [61], where $h_{j}$ accounts for possible molecular
combinations. The stochastic rate constant $c_{j}$ is related to the classical
rate constant $k_{j}$ by $c_{j}=k_{j}V^{1-\nu_{i}}$, where $\nu_{i}$ is the
stoichiometric coefficient. Following Gillespie’s formalism [25], we describe
the chemical Langevin equation (CLE) as follows: If
$\Lambda_{j}[\mathcal{S}(t),\Delta t]$ describes the number of times the
$j^{\mathrm{th}}$ state change occurs in the time interval $(t,t+\Delta t)$,
for $\Delta t>0$ and $\Delta t\ll 1$, then

$\mathcal{S}_{i}(t+\Delta
t)=\mathcal{S}_{i}(t)+\sum_{j=1}^{M}\Lambda_{j}[\mathcal{S}(t),\Delta t]\
\nu_{ji},$ (2)

where $i=1,2,3$ (corresponding to the three molecular species $X$, $Y$, and
$Z$), and $\nu_{ji}$ is a state change vector with $j=1,2,\dots,M$. We now
make important assumptions in Eq. (2). Firstly, in the limit of $\Delta
t\rightarrow 0$, the change in the propensity function $\Delta a\approx 0$.
This holds when the reactant molecular population is large, allowing $\Lambda$
to be approximated by a Poisson random variable, i.e.,
$\Lambda_{j}\rightarrow\mathcal{P}_{j}(a_{j},\Delta t)$. Secondly, in the
limit of $\Delta t\gg 1$ and $a_{j}\gg 1$, one can approximate
$\mathcal{P}_{j}(a_{j},\Delta t)$ by a normal random variable such that
$\mathcal{P}_{j}(a_{j},\Delta t)\rightarrow\mathcal{N}_{j}(a_{j}\Delta
t,a_{j}\Delta t)$. These assumptions are applied simultaneously in the limit
of a large population. With $\mu$ and $\sigma$ as the mean and standard
deviation, respectively, we use
$\mathcal{N}(\mu,\sigma^{2})=\mu+\sigma\mathcal{N}(0,1)$ and put
$\frac{\mathcal{S}_{i}(t+\Delta t)-\mathcal{S}_{i}(t)}{\Delta
t}\approx\frac{d\mathcal{S}_{i}(t)}{dt}$ to arrive at the chemical Langevin
equation (CLE):

$\frac{d\mathcal{S}_{i}(t)}{dt}=\sum_{j=1}^{M}\nu_{ji}\
a_{j}[\mathcal{S}]+\sum_{j=1}^{M}\nu_{ji}\ a_{j}^{1/2}[\mathcal{S}]\ \xi_{j},$
(3)

where the random noise parameters $\xi_{j}=\displaystyle\lim_{dt\rightarrow
0}\mathcal{N}(0,1/dt)$ [25].

###  III.2 Permutation Entropy (PE)

The permutation entropy (PE) method, developed by Bandt and Pompe [58], serves
as a natural complexity measure for time series data. The Bandt-Pompe approach
involves symbolizing a time series by converting it into symbolic sequences.
It determines a probability distribution known as the ordinal probability
distribution for ordinal or permutation patterns in the sequence. The ordinal
probability distribution captures the likelihood of observing different
ordinal patterns, providing valuable insights into the complexity of the time
series. PE is then computed as the Shannon entropy of this ordinal probability
distribution. PE is based on measuring the information contained in comparing
some $r$ consecutive values (known as the order of permutation or embedding
dimension) in a time series.

We explain the Bandt-Pompe procedure in the following. It begins by dividing a
time series $X$ of length $N$, denoted by
$X_{N}=\\{x_{i}~{};~{}i=1,2,\dots,N\\}$, into overlapping partitions
$n=N-(r-1)\tau$, where $\tau$ represents the embedding delay [62]. For each
data partition $D_{p}=(x_{p},x_{p+\tau},\dots,x_{p+(r-1)\tau})$ with
$p=1,\dots,n$ as the partition index, a permutation
$\pi_{p}=(s_{0},s_{1},\dots,s_{r-1})$ of $(0,1,\dots,r-1)$ is determined by
sorting the elements in ascending order. The permutation of the index numbers
is defined by the inequality $x_{p+s_{0}}\leq x_{p+s_{1}}\leq\dots\leq
x_{p+s_{r-1}}$ [58]. The final symbolic sequence is given by
$\\{\pi_{p}\\}_{p=1,\dots,n}$. The relative frequency of all possible patterns
in the symbol sequences is [58]

$\rho_{j}(\pi_{j})=\frac{\hbox{ \\# patterns of type}~{}\pi_{j}~{}\hbox{in
permutation}~{}\\{\pi_{p}\\}}{n}.$ (4)

Then the permutation entropy $S[P]$ is defined as [58]

$S[P]=-\sum_{j=1}^{r!}\rho_{j}(\pi_{j})\log\rho_{j}(\pi_{j}),$ (5)

where the ordinal probability distribution $P=\\{\rho_{j}(\pi_{j})\\}\ ;\
j=1,\dots,r!$. The permutation entropy per symbol of order $r$, denoted by
$h_{r}$, is then given by

$h_{r}=\frac{S[P]}{r-1}.$ (6)

For example, consider the time series $x=(6,2,10,8)$ with an embedding
dimension of $r=3$. Hence, $r!=3!=6$ possible permutations $\\{\pi_{p}\\}$
arise, namely
$\pi_{1}=(0,1,2),\pi_{2}=(0,2,1),\pi_{3}=(1,0,2),\pi_{4}=(1,2,0),\pi_{5}=(2,0,1)$
and $\pi_{6}=(2,1,0)$. Suppose we choose consecutive time units, i.e.,
$\tau=1$. The first data partition $D_{1}=(6,2,10)$ corresponds to
$(x_{t},x_{t+1},x_{t+2})$. Sorting the elements in ascending order yields
$2<6<10$, indicating $x_{t+1}<x_{t}<x_{t+2}$. Therefore, the ordinal pattern
associated with $D_{1}$ is $\pi_{3}=(1,0,2)$. Moving to the second data
partition $D_{2}=(2,10,8)$, sorting the elements in ascending order gives
$x_{t}<x_{t+2}<x_{t+1}$, implying the ordinal pattern is $\pi_{2}=(0,2,1)$.

###  III.3 Complexity-Entropy (CH) causality plane

L$\mathrm{\acute{o}}$pez-Ruiz et al. [59] have introduced the concept of
complexity ($C$) of a system as the product of disequilibrium ($DE$) and
entropy ($H$):

$C=DE\times H.$ (7)

Here, $H$ measures the information the system stores, while $DE$ represents
the system’s deviation from an equiprobable distribution. In essence,
complexity $C$ reflects the interplay between the information stored in a
system (quantified by $H$) and the system’s departure from equiprobability
(quantified by $DE$). We note that entropy measures disorder, while
disequilibrium measures order [63].

Building upon this concept, Rosso et al. [60] have further defined the
statistical complexity measure $C$ as

$C=DE[P,U]\ H[P],$ (8)

where $P$ and $U$ represent the ordinal and uniform probability distributions,
respectively. The uniform distribution is given by $U=\\{1/r!,\ldots,1/r!\\}$.

Connecting with the permutation entropy $S[P]$ defined in Eq. (5), the
normalized permutation entropy $H$ is defined as [60]

$H=\frac{S[P]}{\log r!},$ (9)

where $0\leq H\leq 1$.

In Eq. (8), the disequilibrium $DE$ is expressed as $DE[P,U]=D_{0}D$, where
$D=S\bigg{[}\frac{(P+U)}{2}\bigg{]}-\frac{S[P]}{2}-\frac{S[U]}{2}$ is the
Jensen-Shannon divergence [64] and the normalization constant
$D_{0}=-\frac{1}{2}\bigg{[}\left(\frac{r!+1}{r!}\right)-2\log(2r!)+\log
r!\bigg{]}$ [65, 66].

The complexity-entropy (CH) causality plane, a two-dimensional causality
plane, visually represents the values of the statistical complexity measure
$C$ and the normalized permutation entropy $H$. This plane was used to
distinguish stochastic and chaotic time series [60]. While the normalized
permutation entropy $H$ only measures disorder ($H=0$ for complete order and
$H=1$ for complete disorder), the statistical complexity measure $C$ of Eq.
(8) quantifies both randomness and the degree of correlational structures
[60], making it a powerful tool for analyzing complexity in time series.
Specifically, $C\simeq 0$ indicates both regular and completely random time
series.

##  IV Results and Discussion

We now present the results and discuss our analyses in the following
subsections.

###  IV.1 Regulation of the patterns of cytosolic Ca2+ dynamics by intrinsic
fluctuations

Figure 2: Time evolution of the concentration $x(t)$ of the cytosolic Ca2+ by
solving the chemical Langevin equation (CLE) (12) for the system size: (a)
$V=10^{6}$, (b) $V=10^{4}$, and (c) $V=10^{2}$. Fluctuation-driven dynamics
are observed in the distinct dynamic states of the cytosolic Ca2+, namely,
stable steady-state, simple periodic oscillation, complex oscillations such as
bursting, period doubling, and quasi-periodicity, and chaos. Figure 3: Phase
space plots of the different dynamics of the intracellular Ca2+ oscillation
model (1), namely, (a) simple periodic oscillation, (b) bursting, (c) chaos,
(d) period doubling, and (e) quasi-periodicity for systems of $V=10^{4}$ (blue
curves) and $10^{6}$ (red curves).

To capture the experimentally observed stochastic dynamics of intracellular
Ca2+ oscillations, we adopt a stochastic approach to the non-linear model (1)
using the chemical Langevin equation (3). The stochastic approach will allow
us to analyze the impact of intrinsic fluctuations on the dynamics of
cytosolic Ca2+. Representing $x,y$, and $z$ as the concentrations of cytosolic
Ca2+, stored Ca2+ in the internal pool, and InsP3, respectively, we denote the
state vector of concentrations as
$s=s(t)=[x(t),~{}y(t),~{}z(t)]^{\mathrm{T}}$. We then rewrite the coupled,
non-linear ODEs (1) as

$\frac{ds}{dt}=F(x,y,z),$ (10)

where

$F(x,y,z)=\begin{bmatrix}V_{0}+V_{1}\beta-V_{2}+V_{3}+k_{f}y-kx\\\
V_{2}-V_{3}-k_{f}y\\\ \beta V_{4}-V_{5}-\epsilon z\end{bmatrix}.$ (11)

Using stochastic modeling, we now represent the system of ODEs. (10) as a set
of reaction channels, illustrating the changes of the populations $X,Y$, and
$Z$ (see Table 1 in Appendix A). The state transitions reflect random births
and deaths of the molecular species $X,Y$, and $Z$, introducing intrinsic
fluctuations in their populations. We determine the propensity functions of
each reaction channel using the definition in Sec. III.1. Using Eq. (3), we
determine the chemical Langevin equation (CLE) for the intracellular Ca2+
oscillation model (10) as:

$\frac{ds}{dt}=F(x,y,z)+\frac{1}{\sqrt{V}}\ G(x,y,z),$ (12)

where $V$ denotes the system size, and $G(x,y,z)$ is:

$\displaystyle G(x,y,z)$
$\displaystyle=\begin{bmatrix}\sqrt{V_{0}}\xi_{1}+\sqrt{V_{1}\beta}\xi_{2}-\sqrt{V_{2}}\xi_{3}+\sqrt{V_{3}}\xi_{4}+\sqrt{k_{f}y}\xi_{5}-\sqrt{kx}\xi_{6}\\\
\sqrt{V_{2}}\xi_{7}-\sqrt{V_{3}}\xi_{8}-\sqrt{k_{f}y}\xi_{9}\\\
\sqrt{V_{4}\beta}\xi_{10}-\sqrt{V_{5}}\xi_{11}-\sqrt{\epsilon
z}\xi_{12}\end{bmatrix}.$ (13)

Here, $\xi_{j}$ ($j=1,2,\dots,12$) represents statistically independent
Gaussian white noise with the properties of $\langle\xi_{j}\rangle=0$ and
$\langle\xi_{j}(t)\xi_{j^{\prime}}(t^{\prime})\rangle=\delta_{jj^{\prime}}\delta(t-t^{\prime}).$
While the first term $F(x,y,z)$ in Eq. (12) corresponds to the deterministic
part, the second term represents the stochastic component that introduces
internal fluctuations to the deterministic part. The term $\frac{1}{\sqrt{V}}$
accounts for the effects of system size or intrinsic fluctuation on the Ca2+
oscillation dynamics [67].

We solve the CLE (12) using the Euler method with a total of $10^{4}$
observations and a time step of $0.001$ minutes. Time is measured in minutes
(min), and the concentrations ($x,~{}y,~{}z$) in $\mu$M. Depending on the
parameter values (see Table 2 of Appendix B), cytosolic Ca2+ concentration
$x(t)$ shows various dynamical states, including the stable steady state,
simple periodic oscillation, and complex oscillations such as bursting, period
doubling, quasi-periodicity, and chaos. We refer the readers to Ref. [34] for
detailed characterization of these dynamic states using bifurcation diagrams,
Lyapunov exponents, first-return maps, and power spectra.

Fig. 2 presents the time evolution of the cytosolic Ca2+ concentration, $x(t)$
for three different system sizes: $V=10^{6}$, $10^{4}$, and $10^{2}$.
Correspondingly, the strength of intrinsic fluctuation, measured by
$~{}\frac{1}{\sqrt{V}}$, varies approximately as $\sim 0.001$, $0.01$ and
$0.1$ for $V=10^{6}$, $10^{4}$, and $10^{2}$, respectively. The concentration
$x(t)$ in each panel of Fig. 2 is generated by merging the time series
obtained by solving the CLE (12) under different parameter settings
corresponding to each dynamic state. For a more extended simulation, see Fig.
A.1 in Appendix A. Intrinsic fluctuation regulates the dynamical behavior of
cytosolic Ca2+ concentration depending on the system size.

In Fig. 2(a), the cytosolic Ca2+ concentration $x(t)$ exhibits small
fluctuations, resembling nearly deterministic behavior at system size
$V=10^{6}$. As $V$ becomes large, the CLE (12) approaches
$\frac{ds}{dt}\approx F(x,y,z)$, showing a transition towards the
deterministic limit. Such a transition from stochastic to deterministic
dynamics at large system size is often called the thermodynamic limit [68].

In Fig. 2(b), discernible fluctuations start appearing in the dynamics of
$x(t)$ at $V=10^{4}$. As $V$ decreases, the growing intrinsic fluctuation
significantly impacts the patterns of periodic and chaotic states of $x(t)$.
The periodic doubling state appears to lose its double periodicity and
exhibits behavior akin to chaos at a large value of intrinsic fluctuation.
Previously, intrinsic noise was found to impart chaoticity to periodic limit
cycles in the chemical Lorentz system [69].

In Fig. 2(c), intrinsic fluctuations have a pronounced effect at the small
system size $V=10^{2}$, obliterating the discernible dynamic patterns in
$x(t)$. The large intrinsic fluctuations have obscured the complex
oscillations such as bursting, period doubling, and quasi-periodicity that
were previously differentiable. Chaos becomes no longer distinguishable from
noise.

In Fig. 3, we plot the phase diagrams for various oscillatory dynamics
introduced in Fig. 2 at $V=10^{4}$ (blue) and $10^{6}$ (red). The phase
diagrams reveal distinct cyclic patterns inherent to each dynamic state when
the system is large with negligible intrinsic fluctuations ($V=10^{6}$). We
observe a limit cycle for simple periodic oscillation, a limit cycle with
multiple small loops for bursting, a strange attractor for chaos, a double
loop for period doubling, and a torus for quasi-periodicity. Refer to Fig. A.2
to see chaos and period-doubling clearly. When intrinsic fluctuations become
larger for $V=10^{4}$ (blue), the fine structure in the phase diagram becomes
blurred while the large-scale structure is maintained. If the system size is
further decreased to $V=10^{2}$, noise dominates the deterministic dynamics,
and the phase diagram becomes completely noisy (not shown).

###  IV.2 Interplay of the complexity of cytosolic Ca2+ with intrinsic
fluctuation using permutation entropy

Figure 4: Time evolution of the permutation entropy per symbol of order $r$,
denoted as $h_{r}$, for various dynamic states of cytosolic Ca2+: (a)
$V=10^{6}$, (b) $V=10^{4}$, and (c) $V=10^{2}$. Each panel displays $h_{r}$
calculated for the embedding dimension of $r=3$ (blue), $6$ (orange), and $7$
(green). Vertical red lines mark the distinct dynamic states.

We now explore the complexities in cytosolic Ca2+ concentration using
permutation entropy with the embedding dimension $r$, denoted as $h_{r}$ (see
the definition in Eq. (6)). We consider three values of $r=3$, $6$, and $7$ in
line with practical recommendations [58]. The numerical code for calculating
$h_{r}$ is implemented in Fortran90. The time-evolution of $h_{r}$ is
evaluated for distinct dynamic states of cytosolic Ca2+, i.e., the stable
steady state, simple periodic oscillations, chaos, and complex oscillations
such as bursting, period doubling, and quasi-periodicity. For each dynamic
state, trajectories $x(t)$ are simulated, each with a length of $10^{4}$ and a
step size of $0.001$. The results are presented in Fig. 4, where $h_{r}$
effectively captures the complexities of different dynamic states of cytosolic
Ca2+ concentration and their interplay with intrinsic fluctuations.

In Fig. 4(a), we present the $h_{r}\ (r=3,~{}6,~{}7)$ calculated from the time
series $x(t)$ of cytosolic Ca2+ concentration shown in Fig. 2(a) for
$V=10^{6}$. The permutation entropy $h_{r}$ captures the distinctive features
of different dynamic states in cytosolic Ca2+, with higher values indicating
higher disorder or complexity. The highest value of $h_{r}$ is observed for
the stable steady state. This is attributed to the fact that the steady state
is dominated by random fluctuations without any ordered patterns (see also
Fig. A.1(a)). Following the steady state, the bursting state exhibits moderate
disorder resulting from the fine bursting structures, including spikes with
reduced amplitudes around the plateau. Next to the steady state and bursting,
the chaos state exhibits the third highest $h_{r}$ value. This suggests that a
higher degree of order or regularity is inherent in the chaotic dynamics when
observed through permutation entropy. Lastly, it is noticed that increasing
the permutation order $r$ tends to enhance the resolution, making the
differences in $h_{r}$ between various dynamic states more pronounced.

Continuing to Fig. 4(b), we evaluate $h_{r}$ of $V=10^{4}$ using $x(t)$ in
Fig. 2(b). Compared to the previous case ($V=10^{6}$), the values of $h_{r}$
across all periodic and chaotic states become elevated. The Ca2+ dynamics
attain increased disorder driven by a larger intrinsic fluctuation $(\sim
0.01)$. Particularly noteworthy is the amplified impact of intrinsic
fluctuations when interacting with chaotic or period-doubling states. This
behavior is understood such that chaos generally arises and vanishes via a
period-doubling cascade, obeying a sequence of bifurcations known as
Feigenbaum sequence [34, 70] in non-linear dynamical systems. The observed
higher impact of internal fluctuations on chaos aligns with the previous
report by Wu and Kapral [71] in the Willamowski-R$\mathrm{\ddot{o}}$ssler
model for deterministic chemical chaos [72].

In Fig. 4(c), we plot $h_{r}$ of $V=10^{2}$ using the trajectory $x(t)$ in
Fig. 2(c). In such a small system with strong intrinsic fluctuation, the
$h_{r}$ values of the different dynamic states of the cytosolic Ca2+
concentration are almost the same, signifying a transition from order to
disorder state. As also previously seen in Fig. 2(c), the large internal
fluctuations destroy the distinctive features in the periodic and chaotic
patterns of the cytosolic Ca2+ concentration. Additionally, we compare the
mean permutation entropy $\langle h_{r}\rangle$ among the various Ca2+ dynamic
states for given $r$ and $V$ investigated (Fig. A.3). Similar trends in
$\langle h_{r}\rangle$ for all Ca2+ dynamic states at each $r$ and $V$ are
observed, consistent with the patterns in Fig. 4.

Figure 5: Dynamic states of the cytosolic Ca2+ in the complexity-entropy (CH)
causality plane. Here, the statistical complexity measure $C$ and the
normalized permutation entropy $H$ were calculated using the open-source
Python module ordpy [62]. Six distinct dynamic states are marked on this plane
for $V=10^{2}$ (stars), $10^{4}$ (circles), and $10^{6}$ (pentagons). The
solid lines represent the theoretical curves of the maximum and minimum values
of statistical complexity with the embedding dimension $r=6$. The unfilled
diamonds mark the $(H,C)$ values for fBm processes at several Hurst exponent
$\mathcal{H}$. The black dotted line is the guideline.

###  IV.3 Interplay of the complexity of cytosolic Ca2+ with intrinsic
fluctuation using complexity-entropy ($CH$) causality plane

Figure 6: (a–d) Log-log plots of the normalized permutation entropy $H$ as a
function of $V$ at $r=3,~{}4,~{}5$, and $6$. The markers show the data for the
five distinct cytosolic Ca2+ states. For $V\sim 10^{4}$–$10^{6}$, the $H$
values are fitted (solid lines) with a power-law function $H(V)=aV^{-b}$ where
$a$ is a constant and $b$ is the power-law exponent. (e–h) Statistical
complexity measure $C$ against $V$ at $r=3,~{}4,~{}5,$ and $6$. The pink faded
region indicates the region including the maximum values of $C$,
$C_{\mathrm{max}}$. The green faded region indicates the saturated values of
$C$ in the large system size limit.

We now investigate the interplay of intrinsic fluctuation and the complexity
of various dynamic states of cytosolic Ca2+ on the complexity-entropy (CH)
causality plane for system size $V=10^{6}$, $10^{4}$ and $10^{2}$.

In Fig. 5, we plot the $(H,C)$ components of the six distinct Ca2+ dynamics on
the CH causality plane (where the embedding dimension is $r=6$). The solid
lines represent the theoretical curves for the maximum and minimum of
statistical complexity $C$ [73] (see Appendix B). The analysis shows the six
dynamical states of cytosolic Ca2+ as distinct states on the CH causality
plane, and there exists a nontrivial relationship between statistical
complexity and the entropy of Ca2+ dynamics.

When $V=10^{6}$, all the oscillatory states are placed in the range of
$0<H<0.5$ (see also Fig. B.1), suggesting that the cytosolic Ca2+ dynamics are
in ordered states. As $V$ decreases to $10^{4}$, the ($H,C$) components of all
the periodic and chaotic Ca2+ states shift towards the center of the CH
causality plane, with higher values of $C$ and $H$. This indicates that
compared to $V=10^{6}$, the Ca2+ oscillation system attains more complexity
and information (entropy) driven by a larger intrinsic fluctuation [59]. Upon
further increasing intrinsic fluctuations ($V=10^{2}$), the dynamics of
cytosolic Ca2+ tend to demonstrate larger entropy ($H$) but reduced complexity
($C$). Despite the increase in entropy ($H$), the complexity $C$ decreases
with an increased intrinsic fluctuation due to a reduction in disequilibrium
$DE$, as defined in Eq. (8). In the case of $V=10^{2}$, intrinsic fluctuations
destroy the ordered patterns in all the oscillatory states, signifying a
transition towards disorder states. Importantly, the $V$-dependency shows that
the statistical complexity of cytosolic Ca2+ dynamics attains the maximum when
the system has the intermediate size of $V\sim 10^{4}$.

In Fig. 5, additionally, we compare the complexity-entropy behavior of
cytosolic Ca2+ with that of fractional Brownian motion (fBm) [74, 60]. The
$(H,C)$ component of fBm is evaluated with the Hurst exponent from
$\mathcal{H}=0.5$ to 0.9 [60]. As the comparison shows, intriguingly, the
cytosolic Ca2+ with $V=10^{2}$–$10^{4}$ exhibits a complexity-entropy
structure akin to that of fBm. The $(H,C)$ states of Ca2+ dynamics for
$V=10^{4}$ (circles) follow fBm process with $\mathcal{H}=0.9$ (a persistent
Gaussian random walk [74, 60]). For $V=10^{2}$, notably, the $(H,C)$ values of
Ca2+ dynamics closely align with that of fBm with $\mathcal{H}=0.5$, i.e.,
ordinary Brownian motion.

We further investigate the $H$ and $C$ profiles for varying system size $V$.
Figs. 6(a)–(d) show $H$ vs $V$ in log-log scale for given embedding dimensions
of $r=3$, $4$, $5$, and $6$. We note that $H$ exhibits a power-law behavior in
the intermediate range of $V\sim 10^{4}$–$10^{6}$. The solid lines depict the
best-fit power-law with $H(V)=aV^{-b}$ ($b$: the power-law exponent, $a$: a
constant). The fits have the goodness of $r^{2}$ of $\approx 0.98$–$0.99$, and
the corresponding values of $b$ are given in the legend. This power-law
behavior can be interpreted in such a way that the cytosolic Ca2+ dynamics of
intermediate system size (showing the maximum complexity $C$) contain self-
similar patterns or scale-free dynamical structures. Visual inspection of the
time series $x(t)$ in Figs. 2(a)&(b) signifies such self-similar patterns for
$V=10^{4}$–$10^{6}$. If $V>10^{6}$, as observed, the system under negligible
intrinsic fluctuations exhibits deterministic dynamics without disorder.
Consequently, the measure of disorder, $H$, becomes a constant independent of
$V$ in this regime.

To corroborate the self-similar pattern in the cytosolic Ca2+ dynamics, we
analyze their multifractal structure using the method of multifractal
detrended fluctuation analysis (MFDFA) [75, 76]. This methodology generalizes
the detrended fluctuation analysis (DFA) method [77], numerically identifying
long-range correlations in non-stationary time series. Further description of
the MFDFA algorithm is presented in Appendix C. Here is the main ingredient of
the analysis. (1) For each dynamical state of cytosolic Ca2+, we calculate the
fluctuation function $F_{q}(s)$ [Eq. (17)] for varying scale length $s$ with a
given moment $q$ ranging from $-10$ to $10$. (2) We define the generalized
Hurst exponents $h(q)$ in the scaling relation of $F_{q}(s)\sim s^{h(q)}$ and
estimate $h(q)$ from the fitting of $F_{q}(s)$. Fig. C.1(a) shows $h(q)$ vs
$q$ for $V=10^{4}$ (Left), $10^{5}$ (Middle), and $10^{6}$ (Right). Different
colors represent different dynamic states, as indicated in the legend of the
top left panel. The legend holds for all the plots in Fig. C.1. Notably, it is
observed that $h(q)$ has a significant dependence on $q$, illustrating self-
similarity or, more precisely, multifractality. The latter signifies the
different nature of correlations in large and small fluctuations within the
time series $x(t)$ [75, 76]. Additionally, we calculate the classical
multifractal scaling exponent $\tau(q)$ using Eq. (20) and plot it against $q$
(Fig. C.1(b)). The bi-linear nature of the $\tau(q)$ curves is observed, which
indicates intrinsic similarity (multifractality) arising from non-linear
correlations [78].

We also examine whether $h(q)$ for the various states of cytosolic Ca2+ is
consistent with that of fBm. Fig. C.2 presents the comparison between (a) Ca2+
dynamics of $V=10^{4}$ and fBm of $\mathcal{H}=0.9$ (black dashed line) and
between (b) Ca2+ states of $V=10^{2}$ and fBm of $\mathcal{H}=0.5$. Since fBm
is a monofractal process, as expected, $h(q)$ is evaluated to be a
$q$-independent constant from our MFDFA analysis. Although the monofractal fBm
cannot explain $h(q)$ of the Ca2+ dynamics over the entire $q$ domain, it
explains well the $h(q)$ with $q<0$, i.e., the short-time structures of the
Ca2+ dynamics for both cases of $V=10^{2}$ and $10^{4}$. This finding
implicates that the small-scale or local structure of the Ca2+ time series
resembles the fBm of a specific Hurst index, resulting in the similarity in
the CH causality plane. Note that the large-scale structure of the Ca2+
dynamics, characterized by $h(q>0)$, evidently deviates from that of fBm, as
the cytosolic Ca2+ oscillates on a large-time scale.

In Fig. 6(e)–(h), we plot the statistical complexity $C$ against $V$. For the
embedding dimension $r=3$ [Fig. 6(e)], the statistical complexities
monotonically increase with $V$ and then become stationary. In the cases of
larger embedding dimensions ($r=4$, 5, and 6), the statistical complexities
have the maximum value (indicated by the shaded pink region) at an
intermediate system size $V\sim 10^{4}$–$10^{6}$ and then decrease to reach
stationary values (indicated by the shaded green region). Note that the
intermediate-sized system of $V\sim 10^{4}$–$10^{6}$ in which the statistical
complexity reaches the maximum also exhibits the self-similar pattern in the
cytosolic Ca2+, as we studied above.

Finally, we examine the effect of the embedding dimension $r$ on the
statistical complexity measure $C_{\mathrm{max}}$. Fig. 7 shows
$C_{\mathrm{max}}$ for all the oscillatory cytosolic Ca2+ states for
increasing $r$ up to the largest value that can be investigated within the
total observation time of simulated time series. We find that, in the explored
range of $r$, $C_{\mathrm{max}}$ occurs within the intermediate system size
and monotonically increases with $r$. The latter behavior emerges due to the
fact that the two components of $C$, i.e., permutation entropy and
disequilibrium, increase with the embedding dimension (or permutation order)
$r$, as defined in Eq. (8).

Figure 7: Dependence of the maximum value of the statistical complexity
measure, $C_{\mathrm{max}}$ on the embedding dimension $r$.

##  V concluding remark

The collective behavior of a network (e.g., chemical reaction network) of
diverse components (e.g., molecular species) within a complex system often
leads to a variety of complex dynamic states [79, 80], each playing a
significant functional role. In dynamical systems, a stable steady state
typically signifies an equilibrium state [2], oscillations reflect the active
state of the system [44], and chaos is known to play a crucial role in
information processing [81, 82]. Experimental observations indicate that
calcium ions (Ca2+) within biological cells exhibit a variety of complex
dynamic behaviors, including complex oscillations. The Ca2+ oscillations show
stochastic dynamics and a critical determinant factor for the dynamics is the
system size $V$ of the biological cell. $V$ modulates the intrinsic
fluctuations in Ca2+ dynamics such that intrinsic fluctuation scales as
$\sim\frac{1}{\sqrt{V}}$. A systematic exploration of the interplay between
intrinsic fluctuation and the complexity of observed dynamics of intracellular
Ca2+ has not been thoroughly addressed. Our study addresses this gap through a
comprehensive analysis of the complexities associated with the different
dynamics of cytosolic Ca2+ and their interplay with intrinsic fluctuations.
Leveraging complexity measures such as permutation entropy and statistical
complexity, our results unravel intricate relationships between intrinsic
fluctuations and the complexity of Ca2+ dynamics across various states.

We perform stochastic simulations employing the chemical Langevin equation
(CLE) to model the complex dynamics of intracellular Ca2+. By varying the
system size $V$, we delve into the interplay between intrinsic fluctuation and
the complexity of cytosolic Ca2+ dynamics. Strong intrinsic fluctuations lead
to the breakdown of the ordered states in cytosolic Ca2+ concentration. We
find that permutation entropy effectively characterizes the complexities of
the different dynamic states and their changes due to intrinsic fluctuation.
With permutation entropy, chaos is found to be highly sensitive to intrinsic
fluctuation. The complexity-entropy (CH) causality plane, initially proposed
for distinguishing noise and chaos, proved valuable for assessing the
complexities of diverse, dynamic states and analyzing their interplay with
intrinsic fluctuation. The distinct dynamic states of cytosolic Ca2+ exhibit
varying positions within the theoretical bounds of the CH causality plane as
intrinsic fluctuation varies, indicating varying degrees of complexity and
information content. In an intermediate range of $V$ (or intrinsic
fluctuation), the normalized permutation entropy $H$ follows a power-law
behavior, suggesting the presence of scale-free or self-similar patterns in
the dynamical structures of cytosolic Ca2+, furthermore corroborated by the
multifractal detrended fluctuation analysis (MFDFA), providing additional
insights into the nature of correlations of fluctuations in the time series of
cytosolic Ca2+ concentration. Additionally, we observe peak values in the
statistical complexities of the periodic and chaotic states of cytosolic Ca2+
at an intermediate level of intrinsic fluctuation when adjusting the
permutation order or embedding dimension $r$. This intermediate range is
characterized by self-similar patterns with permutation entropy analysis. Such
high-complexity states may correspond to optimal Ca2+ dynamics, holding
potential biological significance, for instance, information transfer within
signaling pathways. This study deepens our understanding of how intrinsic
fluctuations dynamically regulate the complex behaviors of intracellular Ca2+
across diverse states. Our study is intimately related to understanding the
influence of intrinsic fluctuation on the dynamics of several biological
systems at microscopic and mesoscopic scales. Investigating the effect of
intrinsic fluctuations on the complex dynamics of non-linear biochemical
systems remains a topic of fundamental importance.

Oscillations are well-established far-from-equilibrium phenomena [15]. For a
non-equilibrium system with oscillatory dynamics, the stochastic process along
the trajectory of the oscillating state can be considered as a Brownian motion
along the trajectory, which is dynamically irreversible [83] when
transitioning from an initial state $x_{i}$ to some final state $x_{f}$ and
then from $x_{f}$ to $x_{i}$, or when moving along any two realizations of the
trajectory with the same initial condition. The energy dissipation in such a
non-equilibrium system can be quantitatively measured by estimating the total
entropy production [84], denoted as $\Delta\Sigma_{\mathrm{tot}}$, in the
trajectory. $\Delta\Sigma_{\mathrm{tot}}$ can be calculated from the
probabilities of the forward and backward paths, given by $p[x(t)|x_{i}]$ and
$p[x(t)|x_{f}]$, respectively, such that
$\Delta\Sigma_{\mathrm{tot}}=\ln\frac{p[x(t)|x_{i}]}{p[x(t)|x_{f}]}+\ln\frac{p_{0}(x_{i})}{p_{0}(x_{f})}$,
with normalized distributions $p_{0}(x_{i})$ and $p_{0}(x_{f})$ [85, 86].
Estimating entropy production in stochastic trajectories from experiments or
simulations is an active area of research bearing fundamental importance for a
deeper understanding of non-equilibrium fluctuations and non-equilibrium
properties of dynamical systems. It will be interesting to analyze entropy
production in the non-equilibrium system of intracellular calcium oscillations
and will be carried out in future work.

## Acknowledgements

This research is supported by ALC’s appointment to the YST program at the
APCTP through the Science and Technology Promotion Fund and Lottery Fund of
the Korean Government and the National Research Foundation of Korea, Grant No.
RS-2023-00218927 (JHJ). ALC acknowledges Dr. Mohammad Zubbair Malik
(Scientist, Dasman Diabetes Institute, Kuwait) for the initial discussions.
ALC expresses a special thanks to Yeongjin Kim for his kind help in preparing
the figures.

## Conflict of Interests

The authors have no conflicts to disclose.

## Data Availability Statement

The data that support the findings of this study are available within the
article itself.

## Appendix A

Table 1: Ordinary differential equations describing the intracellular calcium oscillation model (1) with the corresponding reaction channels. The concentrations of the cytosolic Ca2+, stored Ca2+ and InsP3 are represented by $x(=X/V)$, $y(=Y/V)$, and $z(=Z/V)$, respectively, with $V$ as the system size. Ordinary differential equations |  Transition of states |  Propensity functions
---|---|---
$\frac{dx}{dt}=V_{0}+V_{1}\beta-V_{2}+V_{3}+k_{f}y-kx$ |  $X\rightarrow X+1$ |  $V_{0}$
|  $X\rightarrow X+1$ |  $V_{1}\beta$
|  $X\rightarrow X-1$ |  $V_{2}$
|  $X\rightarrow X+1$ |  $V_{3}$
|  $X\rightarrow X+1$ |  $k_{f}y$
|  $X\rightarrow X-1$ |  $kx$
$\frac{dy}{dt}=V_{2}-V_{3}-k_{f}y$ |  $Y\rightarrow Y+1$ |  $V_{2}$
|  $Y\rightarrow Y-1$ |  $V_{3}$
|  $Y\rightarrow Y-1$ |  $k_{f}y$
$\frac{dz}{dt}=\beta V_{4}-V_{5}-\epsilon z$ |  $Z\rightarrow Z+1$ |  $\beta V_{4}$
|  $Z\rightarrow Z-1$ |  $V_{5}$
|  $Z\rightarrow Z-1$ |  $\epsilon z$
Table 2: Parameter values used in the numerical simulation of the intracellular calcium oscillation model (1) [34] Parameters |  Steady state |  Simple |  Bursting |  Chaos |  Period doubling |  Quasi-
---|---|---|---|---|---|---
|  |  periodic |  |  |  |  periodicity
|  |  oscillation |  |  |  |
$V_{0}$ ($\mu$M min-1) |  2 |  2 |  $2$ |  2 |  2 |  2
$V_{1}$ ($\mu$M min-1) |  2 |  2 |  $2$ |  2 |  2 |  2
$\beta$ |  0.01 |  0.5 |  $0.46$ |  0.65 |  0.7 |  0.51
$V_{M2}$ ($\mu$M min-1) |  6 |  6 |  $6$ |  6 |  6 |  6
$k_{2}$ ($\mu$M) |  0.1 |  0.1 |  $0.1$ |  0.1 |  0.1 |  0.1
$V_{M3}$($\mu$M min-1) |  20 |  20 |  $20$ |  30 |  30 |  $20$
$k_{x}$ ($\mu$M) |  0.5 |  0.5 |  $0.3$ |  0.6 |  0.6 |  $0.5$
$k_{y}$ ($\mu$M) |  0.2 |  0.2 |  $0.2$ |  0.3 |  0.3 |  $0.2$
$k_{z}$ ($\mu$M) |  0.2 |  $0.2$ |  $0.1$ |  $0.1$ |  0.1 |  $0.2$
$V_{M5}$($\mu$M min-1) |  $30$ |  $5$ |  $30$ |  $50$ |  50 |  $30$
$k_{5}$ ($\mu$M) |  $0.3$ |  $1$ |  $1$ |  $0.3194$ |  0.3194 |  $0.3$
$k_{d}$ ($\mu$M) |  $0.5$ |  $0.4$ |  $0.6$ |  $1$ |  1 |  $0.5$
$k_{f}$ (min-1) |  $1$ |  $1$ |  $1$ |  $1$ |  1 |  $1$
$k$ (min-1) |  $10$ |  $10$ |  $10$ |  $10$ |  10 |  $10$
$\epsilon$ (min-1) |  $0.1$ |  $0.1$ |  $0.1$ |  $13$ |  13 |  $0.1$
$V_{4}$($\mu$M min-1) |  $5$ |  $2$ |  $2.5$ |  $3$ |  3 |  $5$
$m$ |  2 |  2 |  $4$ |  2 |  2 |  $2$
$p$ |  2 |  2 |  $1$ |  1 |  1 |  $2$
$n$ |  4 |  4 |  $2$ |  4 |  4 |  $4$
Figure A.1: The time evolution of the cytosolic Ca2+ concentration $x(t)$ for
the different dynamic states: (a) stable steady state, (b) simple periodic
oscillation, (c) bursting, (d) chaos, (e) period doubling, and (f) quasi-
periodicity, obtained by solving the chemical Langevin equation (CLE) (12) at
$V=10^{5}$.  Figure A.2: Phase space plots of chaos and period-doubling states
of the intracellular Ca2+, obtained by solving the chemical Langevin equation
(CLE) (12) at the system size $V=10^{10}$. Figure A.3: The mean permutation
entropy $\langle h_{r}\rangle$ per symbol of order $r$ for various dynamic
states of cytosolic Ca2+ at (a) $V=10^{6}$, (b) $V=10^{4}$, and (c)
$V=10^{2}$. For each data point on the plot, the average was determined from
600 time series realizations ($h_{r}$) of 10 minutes each after removing the
initial 500 data points to avoid transient effects. The error bars represent
standard errors.

In Table 1, we perform stochastic modeling of the coupled, non-linear ordinary
differential equations (ODEs) of the intracellular calcium oscillation model
(1). We translate the ODEs into a set of chemical reaction channels that show
the transition of states of the $X,Y$, and $Z$ populations. Table 2 presents
the values of the parameters used to solve the chemical Langevin equation
(CLE) (12) of the intracellular calcium oscillation. The values are taken from
Ref. [34].

## Appendix B

Figure B.1: Location of the different dynamic states of the cytosolic Ca2+ on
the complexity-entropy (CH) causality plane at system size $V=10^{6}$.

For a given normalized permutation entropy $H$, there exists a set of
statistical complexity $C$ values between $C_{\mathrm{min}}$ and
$C_{\mathrm{max}}$ [60]. Here we briefly explain how to calculate
$C_{\mathrm{min}}$ and $C_{\mathrm{max}}$ as described in Ref. [73]. Suppose a
system has $M(=r!)$ possible accessible states $\\{x_{i}\ ;\ i=1,2,\dots,M\\}$
at a given scale with the probability $g_{i}$ of being in state $i$. At
equilibrium, all states are equiprobable with $g_{\mathrm{eq}}=1/M$. For a
given $M$, a set of distributions giving $C_{\mathrm{max}}$ can be
$\\{g_{1},g_{i}\\}$ with $g_{1}=g_{\mathrm{max}}$ and
$g_{i}=\frac{1-g_{\mathrm{max}}}{M-1}$. The index $i=2,3,\dots,M$, and
$g_{\mathrm{max}}$ runs from $1/M$ to 1. Similarly, a set of distributions
giving $C_{\mathrm{min}}$ can be $\\{g_{1},g_{i}\\}$ with
$g_{1}=g_{\mathrm{min}}$ and $g_{i}=\frac{1-g_{\mathrm{min}}}{M-k-1}$, where
$g_{\mathrm{min}}$ runs from 0 to $1/(M-k)$ with $k=0,1,\dots,M-2$.

## Appendix C

Figure C.1: Plots of the multifractal parameters, namely (a) generalized Hurst
exponent $h(q)$, and (b) classical multifractal scaling exponent $\tau(q)$ at
the system size $V=10^{4}$ (Left), $V=10^{5}$ (Middle)and $V=10^{6}$ (Right).
The different dynamical states are depicted by distinct colored curves, as
seen in legend. Figure C.2: Comparison of the generalized Hurst exponent
$h(q)$ of the different dynamical states of Ca2+ with that of fractional
Brownian motion (fBm): the different states of Ca2+ (a) at $V=10^{4}$ with fBm
at $\mathcal{H}=0.9$, and (b) at $V=10^{2}$ with fBm at $\mathcal{H}=0.5$. The
distinct colored curves in (a) and (b) represent the different dynamical
states (see legend in Fig. C.1). The black dashed line represents the
$h(q)\approx\mathrm{constant}$ for fBm, indicating monofractality.

Consider a time series $x_{k}$ of length $N$. The multifractal detrended
fluctuation analysis (MFDFA) algorithm and the associated physical
interpretations are described as follows.

  1. 1.

Calculate the profile of the time series as

$Y(i)\equiv\sum_{k=1}^{i}(x_{k}-\langle x\rangle),\ \ i=1,\dots,N.$ (14)

  2. 2.

Divide the profile $Y(i)$ into
$N_{s}\equiv\mathrm{int}\left(\frac{N}{s}\right)$ non-overlapping segments of
equal scale length $s$. If $N$ is not a multiple of $s$, the same procedure is
repeated from the end resulting in $2N_{s}$ segments.

  3. 3.

Calculate the trend of each $2N_{s}$ segment using the least-square fitting
method. The variance is

$F^{2}(\nu,s)\equiv\frac{1}{s}\sum_{i=1}^{s}\\{Y[(\nu-1)s+i]-y_{\nu}(i)\\}^{2},$
(15)

for each segment $\nu=1,...,N_{s}$, and

$F^{2}(\nu,s)\equiv\frac{1}{s}\sum_{i=1}^{s}\\{Y[(N-(\nu-
N_{s})s+i]-y_{\nu}(i)\\}^{2}$ (16)

for $\nu=N_{s}+1,...,2N_{s}$. The $y_{\nu}(i)$ is the fitting polynomial in
segment $\nu$.

  4. 4.

Determine the $q^{\mathrm{th}}$ order fluctuation function $F_{q}(s)$ by
averaging over all the segments as

$F_{q}(s)\equiv\left\\{\frac{1}{2N_{s}}\sum_{\nu=1}^{2Ns}[F^{2}(\nu,s)]^{q/2}\right\\}^{1/q}.$
(17)

Calculate $F_{q}(s)$ for different time scales $s$ and $q$.

  5. 5.

If the time series $x_{k}$ is long-range power-law correlated, $F_{q}(s)$
follows a power-law with the scale length $s$ as,

$F_{q}(s)\sim s^{h(q)},$ (18)

where the power-law exponent $h(q)$ is known as the generalized Hurst
exponent. In the log-log plots of $F_{q}(s)$ versus $s$ for different values
of $q$, the exponent $h(q)$ corresponds to the slopes of the graphs. For
$q=0$, the fluctuation function $F_{q}(s)$ in Eq. (17) diverges and hence, we
use the expression

$F_{0}(s)\equiv\exp\left\\{\frac{1}{4N_{s}}\sum_{\nu=1}^{2N_{s}}\ln[F^{2}(\nu,s)]\right\\}\sim
s^{h(0)}.$ (19)

If $h(q)$ is a constant independent of $q$, the time series has a monofractal
structure. If the characteristics of small- and large-scale fluctuations
differ, there is a significant dependence of $h(q)$ on $q$, indicating
multifractal behavior. For $q>0$, $h(q)$ describes the scaling behavior of
segments with large fluctuations. For $q<0$, $h(q)$ describes the scaling
behavior of segments with small fluctuations. While $q>0$ accounts for large-
scale or global patterns in time series, $q<0$ accounts for small-scale or
local patterns. The generalized exponent $h(q)$ is related to the classical
multifractal scaling exponent $\tau(q)$ from the standard partition function-
based multifractal formalism through the relation,

$\tau(q)=qh(q)-1.$ (20)

If $h(q)$ has a monofractal behavior, $\tau(q)$ is linear. A non-linear
$\tau(q)$ thus indicates multifractal behavior.

## References

  * [1] M. Mitchell, Complexity: A guided tour, Oxford university press (2009).
  * [2] S.H. Strogartz, Nonlinear dynamics and chaos: With applications to physics, biology, Chemistry and Engineering 441 (1994) .
  * [3] G. Nicolis and I. Prigogine, Exploring complexity: An introduction, .
  * [4] G. Nicolis and C. Nicolis, Foundations of complex systems: emergence, information and predicition, World Scientific (2012).
  * [5] A.N. Pisarchik and U. Feudel, Control of multistability, Physics Reports 540 (2014) 167.
  * [6] D. Rickles, P. Hawe and A. Shiell, A simple guide to chaos and complexity, Journal of Epidemiology & Community Health 61 (2007) 933.
  * [7] J. Aguirre, R.L. Viana and M.A. Sanjuán, Fractal structures in nonlinear dynamics, Reviews of Modern Physics 81 (2009) 333.
  * [8] P.C. Bressloff, Stochastic processes in cell biology, vol. 41, Springer (2014).
  * [9] C. Allen and C.F. Stevens, An evaluation of causes for unreliability of synaptic transmission., PNAS 91 (1994) 10380.
  * [10] V. Shahrezaei and P.S. Swain, The stochastic nature of biochemical networks, Current opinion in biotechnology 19 (2008) 369.
  * [11] R. Blossey, L. Cardelli and A. Phillips, A compositional approach to the stochastic dynamics of gene networks, in Transactions on Computational Systems Biology IV, pp. 99–122, Springer (2006).
  * [12] A.B. Kolomeisky, Motor proteins and molecular motors: how to operate machines at the nanoscale, Journal of Physics: Condensed Matter 25 (2013) 463101.
  * [13] H. Qian, Phosphorylation energy hypothesis: open chemical systems and their biological functions, Annu. Rev. Phys. Chem. 58 (2007) 113.
  * [14] H. Qian, Stochastic physics, complex systems and biology, Quantitative Biology 1 (2013) 50.
  * [15] I. Prigogine, Time, structure, and fluctuations, Science 201 (1978) 777.
  * [16] I.R. Epstein and J.A. Pojman, An introduction to nonlinear chemical dynamics: oscillations, waves, patterns, and chaos, Oxford university press (1998).
  * [17] A. Goldbeter, Dissipative structures in biological systems: bistability, oscillations, spatial patterns and waves, Philos. Trans. R. Soc. A 376 (2018) 20170376.
  * [18] B.J. Chung, B. De Bari, J. Dixon, D. Kondepudi, J. Pateras and A. Vaidya, On the thermodynamics of self-organization in dissipative systems: Reflections on the unification of physics and biology, Fluids 7 (2022) 141.
  * [19] L. Poon and C. Grebogi, Controlling complexity, Physical review letters 75 (1995) 4023.
  * [20] A. Hilfinger and J. Paulsson, Separating intrinsic from extrinsic fluctuations in dynamic biological systems, PNAS 108 (2011) 12167.
  * [21] C.V. Rao, D.M. Wolf and A.P. Arkin, Control, exploitation and tolerance of intracellular noise, Nature 420 (2002) 231.
  * [22] M. Kaern, T.C. Elston, W.J. Blake and J.J. Collins, Stochasticity in gene expression: from theories to phenotypes, Nature Reviews Genetics 6 (2005) 451.
  * [23] M.S. Samoilov, G. Price and A.P. Arkin, From fluctuations to phenotypes: the physiology of noise, Science’s STKE 2006 (2006) re17.
  * [24] D.T. Gillespie, Stochastic simulation of chemical kinetics, Annu. Rev. Phys. Chem. 58 (2007) 35.
  * [25] D.T. Gillespie, The chemical langevin equation, The Journal of Chemical Physics 113 (2000) 297.
  * [26] B. Ye, Y. Zhang, W. Song, S.H. Younger, L.Y. Jan and Y.N. Jan, Growing dendrites and axons differ in their reliance on the secretory pathway, Cell 130 (2007) 717.
  * [27] P. Canham and A.C. Burton, Distribution of size and shape in populations of normal human red cells, Circulation research 22 (1968) 405.
  * [28] Y.-c. Shen, A.N. Shami, L. Moritz, H. Larose, G.L. Manske, Q. Ma et al., Tcf21+ mesenchymal cells contribute to testis somatic cell development, homeostasis, and regeneration in mice, Nature communications 12 (2021) 3876.
  * [29] G. Balázsi, A. Van Oudenaarden and J.J. Collins, Cellular decision making and biological noise: from microbes to mammals, Cell 144 (2011) 910.
  * [30] Y. Li, M. Yi and X. Zou, The linear interplay of intrinsic and extrinsic noises ensures a high accuracy of cell fate selection in budding yeast, Scientific reports 4 (2014) 5764.
  * [31] J.-D.J. Han, Understanding biological functions through molecular networks, Cell research 18 (2008) 224.
  * [32] M. Ackermann, A functional perspective on phenotypic heterogeneity in microorganisms, Nature Reviews Microbiology 13 (2015) 497.
  * [33] C.A. Del Negro, C.G. Wilson, R.J. Butera, H. Rigatto and J.C. Smith, Periodicity, mixed-mode oscillations, and quasiperiodicityin a rhythm-generating neural network, Biophysical Journal 82 (2002) 206.
  * [34] G. Houart, G. Dupont and A. Goldbeter, Bursting, chaos and birhythmicity originating from self-modulation of the inositol 1, 4, 5-trisphosphate signal in a model for intracellular ca2+ oscillations, Bulletin of mathematical biology 61 (1999) 507.
  * [35] T. Matsu-ura, T. Michikawa, T. Inoue, A. Miyawaki, M. Yoshida and K. Mikoshiba, Cytosolic inositol 1, 4, 5-trisphosphate dynamics during intracellular calcium oscillations in living cells, The Journal of cell biology 173 (2006) 755.
  * [36] M. Perc, M. Rupnik, M. Gosak and M. Marhl, Prevalence of stochasticity in experimentally observed responses of pancreatic acinar cells to acetylcholine, Chaos: An Interdisciplinary Journal of Nonlinear Science 19 (2009) .
  * [37] N.A. Tamarina, A. Kuznetsov, C.J. Rhodes, V.P. Bindokas and L.H. Philipson, Inositol (1, 4, 5)-trisphosphate dynamics and intracellular calcium oscillations in pancreatic $\beta$-cells, Diabetes 54 (2005) 3073.
  * [38] D. Wu, Y. Jia, L. Yang, Q. Liu and X. Zhan, Phase synchronization and coherence resonance of stochastic calcium oscillations in coupled hepatocytes, Biophysical chemistry 115 (2005) 37.
  * [39] M. Collier, G. Ji, Y.-X. Wang and M. Kotlikoff, Calcium-induced calcium release in smooth muscle: loose coupling between the action potential and calcium release, The Journal of general physiology 115 (2000) 653.
  * [40] F. Meng, W. To, J. Kirkman-Brown, P. Kumar and Y. Gu, Calcium oscillations induced by atp in human umbilical cord smooth muscle cells, Journal of cellular physiology 213 (2007) 79.
  * [41] A. Verkhratsky and A. Shmigol, Calcium-induced calcium release in neurones, Cell calcium 19 (1996) 1.
  * [42] M.J. Berridge and G. Dupont, Spatial and temporal signalling by calcium, Current opinion in cell biology 6 (1994) 267.
  * [43] K. Thurley, A. Skupin, R. Thul and M. Falcke, Fundamental properties of ca2+ signals, Biochimica et Biophysica Acta-General Subjects 1820 (2012) 1185.
  * [44] R.E. Dolmetsch, K. Xu and R.S. Lewis, Calcium oscillations increase the efficiency and specificity of gene expression, Nature 392 (1998) 933.
  * [45] J. Humeau, J.M. Bravo-San Pedro, I. Vitale, L. Nunez, C. Villalobos, G. Kroemer et al., Calcium signaling and cell cycle:progression or death, Cell calcium 70 (2018) 3.
  * [46] M.C. Pinto, F.M. Tonelli, A.L. Vieira, A.H. Kihara, H. Ulrich and R.R. Resende, Studying complex system: calcium oscillations as attractor of cell differentiation, Integrative Biology 8 (2016) 130.
  * [47] L. Xu, F. Zhang, E. Wang and J. Wang, The potential and flux landscape, lyapunov function and non-equilibrium thermodynamics for dynamic systems and networks with an application to signal-induced ca2+ oscillation, Nonlinearity 26 (2012) R69.
  * [48] H. Puebla, Controlling intracellular calcium oscillations and waves, Journal of Biological Systems 13 (2005) 173.
  * [49] J.M. Borghans, G. Dupont and A. Goldbeter, Complex intracellular calcium oscillations a theoretical exploration of possible mechanisms, Biophysical chemistry 66 (1997) 25.
  * [50] D.A. Terrar, Calcium signaling in the heart, in Calcium Signaling, M.S. Islam, ed., (Cham), pp. 395–443, Springer International Publishing (2020), DOI.
  * [51] M. Perc, M. Gosak and M. Marhl, Periodic calcium waves in coupled cells induced by internal noise, Chemical Physics Letters 437 (2007) 143.
  * [52] F. Folz, K. Mehlhorn and G. Morigi, Interplay of periodic dynamics and noise: Insights from a simple adaptive system, Physical Review E 104 (2021) 054215.
  * [53] R. Friedrich, J. Peinke, M. Sahimi and M.R.R. Tabar, Approaching complexity by stochastic methods: From biological systems to turbulence, Physics Reports 506 (2011) 87.
  * [54] A. Roli, A. Ligot and M. Birattari, Complexity measures: open questions and novel opportunities in the automatic design and analysis of robot swarms, Frontiers in Robotics and AI 6 (2019) 130.
  * [55] S. Lloyd, Measures of complexity: a nonexhaustive list, IEEE Control Systems Magazine 21 (2001) 7.
  * [56] P. Grassberger, Information and complexity measures in dynamical systems, in Information dynamics, pp. 15–33, Springer (1991).
  * [57] C.E. Shannon, A mathematical theory of communication, The Bell system technical journal 27 (1948) 379.
  * [58] C. Bandt and B. Pompe, Permutation entropy: a natural complexity measure for time series, Physical review letters 88 (2002) 174102.
  * [59] R. Lopez-Ruiz, H.L. Mancini and X. Calbet, A statistical measure of complexity, Physics letters A 209 (1995) 321.
  * [60] O.A. Rosso, H. Larrondo, M.T. Martin, A. Plastino and M.A. Fuentes, Distinguishing noise from chaos, Physical review letters 99 (2007) 154102.
  * [61] D.T. Gillespie, Exact stochastic simulation of coupled chemical reactions, The journal of physical chemistry 81 (1977) 2340.
  * [62] A.A. Pessa and H.V. Ribeiro, ordpy: A python package for data analysis with permutation entropy and ordinal network methods, Chaos: An Interdisciplinary Journal of Nonlinear Science 31 (2021) .
  * [63] N. Smaal and J.R.C. Piqueira, Complexity measures for maxwell–boltzmann distribution, Complexity 2021 (2021) 1.
  * [64] M. Martin, A. Plastino and O. Rosso, Generalized statistical complexity measures: Geometrical and analytical properties, Physica A: Statistical Mechanics and its Applications 369 (2006) 439.
  * [65] A. Kowalski, M. Martin, A. Plastino and O. Rosso, Entropic non-triviality, the classical limit and geometry-dynamics correlations, International Journal of Modern Physics B 19 (2005) 2273.
  * [66] L. Zunino, D. Pérez, M. Martín, A. Plastino, M. Garavaglia and O. Rosso, Characterization of gaussian self-similar stochastic processes using wavelet-based informational tools, Physical Review E 75 (2007) 021115.
  * [67] C.-l. Zhu, Y. Jia, Q. Liu, L.-j. Yang and X. Zhan, A mesoscopic stochastic mechanism of cytosolic calcium oscillations, Biophysical Chemistry 125 (2007) 201.
  * [68] D.T. Gillespie, The chemical langevin and fokker- planck equations for the reversible isomerization reaction, The Journal of Physical Chemistry A 106 (2002) 5063.
  * [69] U.S. Thounaojam, Stochastic chaos in chemical lorenz system: Interplay of intrinsic noise and nonlinearity, Chaos, Solitons & Fractals 165 (2022) 112763.
  * [70] Y. Pomeau and C. Vidal, Order within chaos: towards a deterministic approach to turbulence, Wiley (1984).
  * [71] X.-G. Wu and R. Kapral, Internal fluctuations and deterministic chemical chaos, PRL 70 (1993) 1940.
  * [72] K.-D. Willamowski and O. Rössler, Irregular oscillations in a realistic abstract quadratic mass action system, Zeitschrift für Naturforschung A 35 (1980) 317.
  * [73] X. Calbet and R. López-Ruiz, Tendency towards maximum complexity in a nonequilibrium isolated system, Physical Review E 63 (2001) 066116.
  * [74] J.-H. Jeon and R. Metzler, Fractional brownian motion and motion governed by the fractional langevin equation in confined geometries, Physical Review E 81 (2010) 021103.
  * [75] J.W. Kantelhardt, E. Koscielny-Bunde, H.H. Rego, S. Havlin and A. Bunde, Detecting long-range correlations with detrended fluctuation analysis, Phys. A: Stat. 295 (2001) 441.
  * [76] 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, Phys. A: Stat. 316 (2002) 87.
  * [77] C.-K. Peng, S.V. Buldyrev, S. Havlin, M. Simons, H.E. Stanley and A.L. Goldberger, Mosaic organization of dna nucleotides, Physical review e 49 (1994) 1685.
  * [78] W.-X. Zhou, D. Sornette and W.-K. Yuan, Inverse statistics and multifractality of exit distances in 3d fully developed turbulence, Physica D: Nonlinear Phenomena 214 (2006) 55.
  * [79] S. Boccaletti, V. Latora, Y. Moreno, M. Chavez and D.-U. Hwang, Complex networks: Structure and dynamics, Physics reports 424 (2006) 175.
  * [80] V. Eskov, V. Eskov, J. Vochmina and T. Gavrilenko, The evolution of the chaotic dynamics of collective modes as a method for the behavioral description of living systems, Moscow university physics bulletin 71 (2016) 143.
  * [81] H. Korn and P. Faure, Is there chaos in the brain? ii. experimental evidence and related models, Comptes rendus biologies 326 (2003) 787.
  * [82] A. Riaz and M. Ali, Chaotic communications, their applications and advantages over traditional methods of communication, in 2008 6th International Symposium on Communication Systems, Networks and Digital Signal Processing, pp. 21–24, IEEE, 2008.
  * [83] D.S. Seara, B.B. Machta and M.P. Murrell, Irreversibility in dynamical phases and transitions, Nature communications 12 (2021) 392.
  * [84] U. Seifert, Fluctuation theorem for birth–death or chemical master equations with time-dependent rates, J. Phys. A 37 (2004) L517.
  * [85] D. Andrieux, P. Gaspard, S. Ciliberto, N. Garnier, S. Joubaud and A. Petrosyan, Entropy production and time asymmetry in nonequilibrium fluctuations, Physical review letters 98 (2007) 150601.
  * [86] T. Xiao, Z. Hou and H. Xin, Stochastic thermodynamics in mesoscopic chemical oscillation systems, The Journal of Physical Chemistry B 113 (2009) 9316.