# Privacy-Assured Outsourcing of Compressed Sensing Reconstruction Service in
Cloud

Ping Wang P. Wang was with the College of Electronics and Information
Engineering, Southwest University, Chongqing 400715, China (e-mail:
bruce_wp@163.com).

###### Abstract

Compressed sensing (CS), breaking the constriction of Shannon-Nyquist sampling
theorem, is a very promising data acquisition technique in the era of
multimedia big data. However, the high complexity of CS reconstruction
algorithm is a big trouble for endusers who are hardly provided with great
computing power. The combination of CS and cloud has the potential of freeing
endusers from the resource constraint by cleverly transforming computational
workload from the local cilent to the cloud platform. As a result, the low-
complexity encoding virtue of CS is fully leveraged in the resource-
constrained sensing devices but its highcomplexity decoding problem is
effectively addressed in cloud. It seems to be perfect but privacy and
security concerns are ignored. In this paper, a secure outsourcing scheme for
CS reconstruction service is proposed. Experimental results and security
analyses demonstrate that the proposed scheme can restrict malicious access,
verify the integrity of the recovered data, and resist brute-force attack,
ciphertext-only attack, and plaintext attack.

###### Index Terms:

compressed sensing, sparse signal reconstruction, cloud, outsourcing, security
and privacy.

##  I Introduction

With the development of multimeadia big data, how to acquire data efficiently
is an important topic all the time. In the traditinal infromation acquisition
system, the encoder usually requires to acquire massive amounts of data
following Shannon-Nyquist sampling theorem and then compress it for efficient
usage of storage and bandwidth resources. Such a sampling-then-compression
framework works in a resource-unfriendly manner and would take a great deal of
complexity on the design of sensing system. Compressed sensing (CS) is a
newly-developing data acquisition and reconstruction technique mixing the
signal sampling and data compression processes by virtue of the sparsity of
signal to be captured. In the CS paradigm, the combination of sampling and
compression could not only save a considerable lot of computing and storage
resources but also cut down the hardware overhead sharply. In the past decade,
CS takes a significant step from theory to practice, which is widely used in
magnetic resonance imaging[1], sparse channel estimation[2], and internet of
thighs (IoT)[3], etc. All of these take advantage of the low-complexity
encoding virtue of CS.

Every coin has two sides, and CS is no exception. In the CS-based
applications, the low-complexity encoding is at the cost of the high-
complexity decoding. Accordingly, CS can be called as a magic technique of
shifting the complexity from the encoder to the decoder. In essence, decoding
the captured data, namely sparse signal reconstruction (SSR), is to solve a
linear programming problem which requires a substantial amount of
computational power. For the resource-limited endusers, it is not an advisable
choice to put such a heavy workload on the terminal devices, such as mobile
phone and some wearable IoT devices. On the one hand, the workload of SSR task
could not usually be well-matched with the computational power of the terminal
devices. On the other hand, performing SSR in the local client is a time-
consuming task. Therefore, how to perform the SSR task effectively has always
been a key challenge of restricting the CS-based applications. With the
development of cloud computing, an alternative approach is to outsourced the
SSR task to a cloud.

Cloud Computing permits convenient on-demand network accesses to a shared pool
of configurable resources (e.g., networks, storage, and computations) that can
be rapidly provided and released with great efficiency and minimal management
effort[]. Delivering the abundant resources to the customers is the main
service of cloud providers. By outsourcing the workloads to the cloud, the
computational power of the endusers is no longer limited by their resource-
constrained devices. Through the above discussions, outsourcing SSR services
to a cloud seems to be feasible and promising. However, there still exist a
number of fundamental and critical challenges, among which security is the top
concern.

With the combination of CS and IoT, the outsourced SSR service often contain
some sensitive information, such as military radar data and personal health
monitoring data. Most traditional cryptosystems, providing a high level of
security guarantee at the cost of a heavy computation workload, is hardly
embedded in the physical layer of the lightweight sensing applications.
Therefore, there exist some practical security challenges in such an
outsourced framework where the high-complexity SSR task is transformed from
the resource-constrained sensing devices to the cloud and the reconstructed
data are sent to the legitimate enduers. At first, how to guarantee the
security of transmitting sensitive data from the sensing devices to the cloud
servers is of great concern in the absence of specialized cryptosystems.
Secondly, the reconstructed data should be secret for the semi-trusted clouds
but available for the legitimate endusers. In addition, from the perspective
of a cloud provider, it plays a pivotal role how to check the legality of
enduers asking for an access to the reconstructed data. Finally, how to avoid
the tampering attack to the reconstructed data but permit content-preserving
data distortions is a big challenge all the time, resulting from that the
reconstructed data are doomed to be different from the original data to some
extent.

In this paper, a privacy-assured cloud-assisted outsourcing scheme is proposed
for the SSR service. Without any additional cryptosystem, two privacy-
preserving primitives, compressed sensing-based encryption and sparse
represent-based encryption, are seamlessly embedded in the process of CS
encoding. The cascade effect of two primitives provides a high level of
transmission secrecy guarantee for the captured data. Sparse represent-based
encryption makes sure that the reconstructed data would not leak any
information except for sparsity. To restricting malicious access, a kind of
access passwords is extracted from some prior encoding knowledge and then
regarded as the certificate of authorized endusers. Besides, CS-based message
authentication codes are generated to figure out whether the receiving data
have been pulluted or tampered.

The rest of this paper is organized as follows. Some related works are briefly
introduced and analyzed in Section 2. Section 3 gives the basic knowledge of
CS, system architecture, threat model, and design goals. In Section 4, several
basic techniques are proposed, including bi-level encryption, integrity
verification, and access control. Section 5 describes our privacy-assured
outsourcing scheme in detail. Experimental results and analysis are given in
Section 6 . A conclusion is drawn in Section 7.

##  II Related works

SSR is essentially a linear programing problem and its solution generally
needs the support of a great computational power. Outsourcing the high-
complexity CS decoding task to the cloud is a newly well-studied flied. Secure
computation outsourcing[4, 5, 6] is to protect the input and ouput privacy in
the process of outsourced computation. As a breakthrough work, the fully
homomorphic encryption (FHE) was proposed to enables directly performing
functions on encrypted data while yielding the same encrypted results as if
the functions were run on plaintext. Aftherwards, some secure outsourcing
schemes based on FHE were proposed[7, 8, 9]. However, they are usually
impractical because of the exceedingly high complexity of FHE operation. Note
that SSR workloads mainly consist of matrix transformation and linear equation
computations. Lei et al. proposed a number of secure outsourcing protocols for
solving matrix inversion[10], matrix multiplication[11], and matrix
determinant computations[12]. Wang et al. presented some privacy-preserving
outsourcing protocols for linear equation[13, 14] and liner programming[15,
16] in cloud. Besides, Wang et al. further proposed a secure outsourcing
scheme for CS-based image reconstruction service[17] but its complexity is
still too high for the resource-constrained devices To achieve secure SSR
outsourcing as possible in a low-complexity manner possible, some newly works
tried to transform SSR problem to the encrypted one in the process of SR[18,
19, 20]. All of those aimed to embed cryptographic operation in the process of
spare representation (SR). Such an ideology occurred in some CS-based privacy-
preserving works[21, 22, 23]. However, SR-based encryption is a weak symmetric
cipher such that anyone holding abundant computing resources can reveal the
sparsity of the plaintext. It remain challenging to provide outsourced data
with high level of security guarantee in a resource-constrained scenario.

##  III Problem formulation

This section introduces the basic knowledge of CS, system architecture, threat
model, and our design goals.

###  III-A CS preliminaries

According to CS theory [24, 25, 26], there are two preconditions of
guaranteeing the successful reconstruction from very few samples. One is that
the signal to be sampled should be sparse or could be sparsely represented in
a certain basis. Fortunately, most natural signals are sparse under an
appropriate transformation domain. The other is that the measurement matrix
should be as incoherent with the basis matrix as possible. CS theory can be
mathematically summarized as follows. An one-dimensionality signal ${\bf
x}\in\mathbb{R}{{}^{N}}$ is sparsely represented as a coefficient vector ${\bf
s}\in\mathbb{R}{{}^{N}}$ containing at most $K$ non-zero elements in a certain
basis ${\bf\Psi}\in{\mathbb{R}^{N\times N}}$ as follows:

$\bf x={\bf\Psi}s,$ (1)

where ${\left\|{\bf s}\right\|_{0}}=K\ll N$ and ${\left\|{\cdot}\right\|_{0}}$
counts the number of the non-zero entries. In the CS paradigm, such a
coefficient vector ${\bf s}$ can be sampled and compressed simultaneously by
the sensing matrix ${\bf A}\in{\mathbb{R}^{M\times N}}$ as follows:

${\bf y}={\bf A\Psi}^{-1}{\bf x}={\bf As},$ (2)

where ${\bf y}\in\mathbb{R}{{}^{M}}$ $(K<M\ll N)$ denotes a measurement vector
(i.e., sampled data) and ${\bf\Phi}={\bf A}{\bf\Phi}^{-1}$ denotes the
measurement matrix. ${\bf A}$ must obey the energy-preserving guarantees [27,
28, 29].

Evidently, the linear systems $(2)$ does not have the unique solution.
Considering the sparsity of this signal, the intuitive way is to solve
${l_{0}}$ norm problem known as greedy iteration with denoising as follows:

${\bf s}^{\prime}=\mathop{\arg\min}\limits_{\bf s}{\left\|\bf
s\right\|_{0}}\quad s.t.\quad{\rm{}}{\left\|{\bf
y-As}\right\|_{2}}\leq\varepsilon,$ (3)

where $\varepsilon$ indicates noise levels. In essence, solving $(3)$ is a NP-
hard problem because it works with exhaustively searching over all column
subsets of $\bf{A}$, which can be relaxed as ${l_{1}}$ norm problem [27] known
as basis pursuit with denoising as follows:

${\bf s}^{\prime}=\mathop{\arg\min}\limits_{\bf s}{\left\|\bf
s\right\|_{1}}\quad s.t.\quad{\left\|{\bf y-As}\right\|_{2}}\leq\varepsilon,$
(4)

where ${\left\|{\cdot}\right\|_{1}}$ counts the sum of absolute value of all
elements. The basis pursuit algorithm based on solving $(4)$ is a linear
programming problem, thus it can relatively accurately reconstruct this
original signal at the cost of high computation.

Except for the sparsity, there is another precondition that ${\bf\Phi}$ should
be as incoherent as possible with ${\bf\Psi}$. Strictly speaking, the
incoherence means that any basis vector in ${\bf\Psi}$ can not be represented
in the basis composed of the rows of ${\bf\Phi}$. To meet such a incoherence
requirement, a mathematical definition of Restricted Isometry Property (RIP)
was given in [28], in which a matrix ${\bf\Phi}$ is said to satisfy the RIP of
order $K$ if there exists a ${\delta_{K}}\in(0,1)$ such that

$(1-{\delta_{K}})\left\|\bf s\right\|_{2}^{2}\leq\left\|{{\bf\Phi}\bf
s}\right\|_{2}^{2}\leq(1+{\delta_{K}})\left\|\bf s\right\|_{2}^{2}$ (5)

holds for all $k$-sparse signal $\bf s$.

###  III-B System architecture and threat models

Figure 1: System architecture.

By virtue of low-complexity encoding, CS has a great potential of being
considered as a low-overhead sensing technique widely used in the physical
layer of IoT. With the combination of cloud and IoT, a promising CS-based
application scenario is shown in Fig. 1. In a CS-based sensor network, the
subnodes capture the original data $\bf x$ in a cost-effective manner and then
the parent nodes outsource the SSR task ${\Gamma}=({\bf A},{\bf y})$ to a
remote cloud. Afther receiving the access request from the enduser, the cloud
manager firstly verify the legitimacy of this enduser and then determine
whether or not to send the reconstructed result. In the local side, the
enduser firstly figure out whether the received data ${\bf s}^{\prime}$ haven
been tampered or polluted and then gets the final data ${\bf x}^{\prime}$ by a
simple computation ${\bf x}^{\prime}={\bf\Psi}{\bf s}^{\prime}$ .

In this paper, we consider a semi-trusted cloud. Specifically, the cloud
provider is supposed to honestly perform SSR task but very curious about the
information contained by computational results. There exist three main threats
in the above system architecture, resulting from two deep-rooted practical
issues. One is how to provide the caputured data with high confidentiality in
the resource-constrained sensing player. The other is how to buid the trust
mechanism in the process of sensor-cloud-user interaction. For each sensor
node, most traditional high-energy cryptosystems are unpractical such that the
captured privacy-sensitive data are attractive to the malicious attackers. In
the absence of the privacy-assured measures, senor and user lose the ownship
of the reconstructed result as a result of outsourcing SSR task to the cloud.
In the process of cloud-user interaction, attackers may disguise themselves as
a legitimate enduser to access to the reconstructed result. In addtion, the
atttackers may tamper the sampling data to make SSR fail or distort original
information.

###  III-C Design Goals

In consideration of the privious system architecture and threat models, we
proposed a privacy-assured outsourcing framework for SSR service in cloud to
meet the following design goals:

  * •

Transmission confidentiality: The confidentiality of the captured data is
cleverly guaranteed without a specialized cryptosystem.

  * •

Privacy-assured outsourcing: The semi-trusted cloud has to honestly perform
the SSR task and has no ability of stealing the privacy contained in the
reconstructed data.

  * •

Access control: The malicious access to the reconstructed data can be
efficiently restricted.

  * •

Integrity verification: The energy-infinite noises can be tolerated but the
tampering attack can be exactly detected.

##  IV Basic techniques

In this section, we propose three key techniques, including bi-level
encryption, integrity verification, and access control.

###  IV-A Bi-level encryption

Most existing high-complexity cryptosystems usually could not be matched with
the resource-constrianed information acquisition system such that how to
guarantee the confidentiality of the captured data are an ubiquitous
challenge. In this paper, we propose a bi-level encryption scheme based on CS,
where the cryptographic features are embedded in both subsampling and sparse
representation processes. Specifically, two key-controlled sensing matrix
${\bf\Phi}_{k}$ and basis matrix ${\bf\Psi}_{k}$ are designed, renewed, and
kept secret.

In [21], Zhang at el. indicated that andomly exchanging the position of rows
and multiplying rows by a non zero number would not alter the SR competence of
a certain sparsifying basis, which can be mathematically summarized in Theorem
1.

###### Theorem 1

Two basis matrices ${\bf\Psi}$, ${\bf\Psi}^{\prime}$ are said to be
equivalent, namely ${\bf x}={\bf\Psi s}={\bf\Psi}^{\prime}{\bf s}^{\prime}$,
${\left\|{\bf s}\right\|_{2}}={\left\|{\bf s}^{\prime}\right\|_{2}}$, if given
any permutation matrix ${\bf P}$ and diagonal matrix ${\bf D}$ there exists
${\bf\Psi}^{\prime}={\bf\Psi}{\bf PD}$.

Although the spare representations are different under different sparsifying
basis, their sparsity would be doomed to be same. According to it, the
sparsity representation process can be viewed as a symmetric encryption
process as follows:

${\bf s}=({{\bf\Psi}}{{\bf P}_{k_{1}}}{{\bf D}_{k_{2}}})^{-1}{\bf
x}={{\bf\Psi}^{-1}_{k}}{\bf x}$ (6)

where ${\bf x}$, ${\bf s}$, and ${{\bf\Psi}_{k}}={{\bf\Psi}}{{\bf
P}_{k_{1}}}{{\bf D}_{k_{2}}}$ are regarded as the plaintext, the ciphertext,
and the secret, respectively, as well as ${\bf P}_{k_{1}}$ and ${\bf
D}_{k_{2}}$ denote the key-controlled random permutation matrix and diagonal
matrix, respectively. Accordingly, the semi-trusted cloud providers could not
get the privacy-sensitive ${\bf x}$ from the reconstructed ${\bf s}$ with the
absence of ${{\bf\Psi}_{k}}$.

Some previous works encrypts and sparsify the original signal simultaneously
and then directly outsource the finally acquired data ${\bf y}={\bf A}{\bf s}$
to the cloud[18, 19, 20]. In their outsourcing framework, the sparsity of the
plaintext is a public information for everybody with strong computational
power. From the perspective of cryptography, their outsourcing frameworkes are
vulnerable and fragile. To fix such a security flaw, we try to embed the
cryptographic features in the subsampling process, namely CS-based encryption.

CS-based encryption is a recently proposed lightweight cipher, in which ${\bf
s}$, ${\bf y}$, and ${\bf A}$ are regarded as the plaintext, the ciphertext,
and the secret key, respectively, as well as the subsampling and the
reconstruction processes as the encryption and the decryption processes,
respectively. The key is to make sure that the sensing matrix ${\bf A}$ is
known by the encoder and the decoder and satisfies the RIP condition.
Considering the security and efficiency, we employ the structurally subsampled
matrix to construct ${\bf A}$ controlled by several parameters and satisfying
RIP condition.

At fist, a block identity matrix ${\bf S}\in{\mathbb{R}^{M\times N}}$ is
constructed as defined in Definition 1, and a diagonal matrix ${\bf
R}_{k_{3}}\in{\mathbb{R}^{N\times N}}$ whose non zero elements are independent
symmetric Bernoulli or Rademacher random variables is generated and renewed by
changing the original conditions or inputs of pseudo random number generator
(PRNG). Then, a structurally subsampled matrix is generated as follows:

${\bf A}_{k}={\bf S}{\bf R}_{k_{3}}{\bf\Psi},$ (7)

where $\bf\Psi$ denotes a unitary matrix (i.e., sparsifying basis).
Essentially, the combination form ${\bf S}{\bf R}_{k_{3}}$ is a row-mixing
matrix. According to Theorem 2 given in [30], such a ${\bf A}_{k}$ has been
proved to be a RIP matrix.

###### Deifinition 1

Given $N$ and its divisor $L$, a block identity matrix is defined as ${\bf
S}=[{{S}_{1}^{T}},{{S}_{2}^{T}},...,{{S}_{M}^{T}}]^{T}$, whose $i$-th row
vector $[{{\bf 0}^{1\times(i-1)L}}{{\bf 1}^{1\times L}}{{\bf
0}^{1\times(N-iL)}}]$ containing $L$ adjacent $1$s is $L$-sparse, $i\in[1,M]$,
$M={N\mathord{\left/{\vphantom{NL}}\right.\kern-1.2pt}L}$.

###### Theorem 2

[30] Let ${\bf A}\in{\mathbb{R}^{H\times N}}$ be a structurally subsampled
matrix, randomly subsampled from the unitary matrix $\bf U$ and the row-mixing
matrix ${\bf S}{\bf R}$ and then normalized column by column, namely, ${\bf
A}=\sqrt{\frac{M}{H}}{\bf S}{\bf R}{\bf U}$, $H\leq M$. For each integer
$a>2$, any $b>1$, and any $\delta_{K}\in({\rm{0}},{\rm{1}})$, there exist two
absolute positive constants $c_{1}$ and $c_{2}$ such that the matrix ${\bf A}$
satisfies the RIP of order $K$ with probability at least
${\rm{1-20max\\{exp(-}}{{\rm{c}}_{2}}\delta_{K}^{2}{\rm{z),}}{{\rm{N}}^{-1}}{\rm{\\}}}$
if

${\rm{H}}\geq{c_{1}}zKN{\mu^{2}}({\bf U}){\log^{3}}N{\log^{2}}K,$ (8)

where ${\mu}({\bf U})$ denotes the coherence coefficient of ${\bf U}$, namely

$\mu({\bf U})=\mathop{\max}\limits_{1\leq i\neq j\leq
N}\frac{{\left|{\left\langle{{\psi_{i}},{\psi_{j}}}\right\rangle}\right|}}{{{{\left\|{{\psi_{i}}}\right\|}_{2}}{{\left\|{{\psi_{j}}}\right\|}_{2}}}}.$
(9)

Unifying SR and subsampling processes, the CS projection can be expressed as :

${\bf y}={\bf S}{\bf R}_{k_{3}}{\bf\Psi}({{\bf\Psi}}{{\bf P}_{k_{1}}}{{\bf
D}_{k_{2}}})^{-1}{\bf x}={\bf A}_{k}{\bf s}.$ (10)

The above format implies that the original signal $\bf x$ is sparsely
represented in the transformation domain ${\bf\Psi}_{k}$ and the sensing
matrix ${\bf A}_{k}={\bf S}{\bf R}_{k}{\bf\Psi}$ is generated from the row-
mixing matrix ${\bf S}{\bf R}_{k}$ and ${\bf\Psi}$ in a structurally
subsampling manner. Although there exist some differences between ${\bf\Psi}$
and ${\bf\Psi}_{k}$, ${\bf A}_{k}$ is still a RIP matrix, which can be
inferred from Theorem 3.

###### Lemma 1

For two equivalent basis matrix ${\bf\Psi}$ and ${\bf\Psi}^{\prime}$, their
coherence coefficients are also equivalent, namely
$\mu({\bf\Psi})=\mu({\bf\Psi}^{\prime})=\mu({\bf\Psi},{\bf\Psi}^{\prime})$.

Proof. The equivalence between ${\bf\Psi}$ and ${\bf\Psi}^{\prime}$ means that
there must be a permutation matrix $\bf P$ and a diagonal matrix $\bf D$,
which satisfy ${\bf\Psi}^{\prime}={\bf\Psi}{\bf PD}$. Aiming at column
permutations, $\bf P$ obviously does not affect the coherence. $\bf D$ aims to
multiply some rows by non zero numbers. It makes the inner product and 2-norm
increase by equal times such that the coherence still retains unchanged. All
in all, ${\bf\Psi}$ and ${\bf\Psi}^{\prime}$ do not the coherence and the
above lemma holds.

###### Theorem 3

If a structurally subsampled matrix ${\bf A}$, generated from an unitary
matrix $\bf U$ and a row-mixing matrix ${\bf S}{\bf R}$, satisfies the RIP
with high probability, a matrix ${\bf A}^{\prime}$ satisfying RIP with equally
high probability can be generated from ${\bf U}^{\prime}={\bf U}{\bf P}{\bf
D}$ and ${\bf S}{\bf R}$, where ${\bf P}$ and ${\bf D}$ denote the permutation
matrix and the diagonal matrix, respectively.

Proof. The above theorem implies that permutation matrix and diagonal matrix
do not affect the RIP condition for any structurally subsampled matrix. $1$)
${\bf A}$ is generated by subsampling ${\bf SR}{\bf U}$ and then normalizing
its column vectors. The normalize operation makes ${\bf D}$ in vain. So
subsampling ${\bf U}{\bf P}{\bf D}$ is equivalent to subsampling ${\bf U}{\bf
P}$. Furthermore, ${\bf U}{\bf P}$ is also an unitary matrix.

According to Theorem 2, ${\bf A}^{\prime}$ satisfy RIP of order $K$ with
probability at least
${\rm{1-20max\\{exp(-}}{{\rm{c}}_{2}}\delta_{K}^{2}{\rm{z),}}{{\rm{N}}^{-1}}{\rm{\\}}}$
if

${\rm{H}}\geq{c_{1}}zKN{\mu^{2}}({\bf U}{\bf P}{\bf
D}){\log^{3}}N{\log^{2}}K.$ (11)

According to Lemma 1, ${\mu}({\bf U}{\bf P}{\bf D})={\mu}({\bf U})$ ensures
that $(11)$ is equivalent to $(9)$. So the above theorem holds.

By keeping the keys $(k_{1},k_{2},k_{3})$ secretly renewed between the encoder
and the decoder, such a projection like $(10)$ is a bi-level symmetric
encryption process. $(k_{1},k_{2})$ are regarded as the secret keys of SR-
based encryption and $k_{3}$ as the secret key of CS-based encryption.

###  IV-B Integrity verification

Outsourced computations involve multi individuals, including sensor, cloud,
and enduser. Due to the presenceof maliciou attackers and the absence of trust
mechanism, data may be polluted or tampered in the process of transmission,
storage, and computing. So guaranteeing data integrity is of great importance.
Unfortunately, CS data integrity verification is a difficult work in the
resource-limited scenarios. One is how to generate message authentication
codes $MAC$ directly for a analog signal. The other is how to distinguish
malicious tamper and acceptable noise in view of the robustness of CS
measurements.

It is well known that CS is a lossy and lightweight data acquisition technique
such that most existing integrity verification algorithms like MD5 and SHA2,
which are very sensitive for data perturbation and endowed with high
complexity, are infeasible in the CS-based applications. Therefore, it is
challenging to verify the integrity of the received data in the CS scenarios.

By leveraging the compressive feature and breaking RIP condition, an one-way
hash matrix ${\bf A}_{MAC}\in{\mathbb{R}^{m\times N}}$ $(m<<N)$ can be
constructed to get the authentication measurements ${\bf y}_{MAC}$ of sparse
coefficients $\bf s$ while subsampling as follows:

${{\bf y}_{MAC}}={{\bf A}_{MAC}}{{\bf\Psi}^{-1}}{\bf x}={{\bf A}_{MAC}}{\bf
s}.$ (12)

Note that the one-way guarantee of the above process depends on the non-RIP
condition of ${\bf A}_{MAC}$. Considering that, we suggest that ${\bf
A}_{MAC}={\bf\Phi}_{MAC}{\bf\Psi}$ is a non-RIP matrix by making
${\bf\Phi}_{MAC}$ as incoherent with ${\bf\Psi}$ as possible and letting $m$
less than what is required by the RIP condition.

As studied in [31], the relation between the incoherence and the number of the
required measurements is given in Theorem 4. Accordingly, it is not hard to
draw a conclusion that the matrix ${\bf A}_{MAC}$ is one-way as long as its
rows number is less than $CK\sqrt{N}t\mu({\Phi^{T}},\Psi)\log(tK\log
N){\log^{2}}K$.

###### Theorem 4

[31] For a fixed basis matrix $\bf\Psi$, the matrix
$\bf\Phi\Psi\in{\mathbb{R}^{M\times N}}$ is said to satisfy the RIP condition
with probability as least $1-5{e^{-t}}$ if

$M\geq CK\sqrt{N}t\mu({\Phi^{T}},\Psi)\log(tK\log N){\log^{2}}K$ (13)

To relieve the transmission burden and increase the entropy, $MAC$ are
generated by extracting the sign of ${{\bf y}_{MAC}}$ and then performing non-
return-to-zero (NRZ) line code as follow:

$MAC={f_{NRZ}}(sign({{\bf{y}}_{MAC}}))={f_{NRZ}}(sign({{\bf{A}}_{MAC}}{\bf{s}}),$
(14)

where ${f_{NRZ}}(\cdot)$ denotes a NRZ encoder. Here, we suppose that the
security of transmitting ${\rm MAC}$ from the sender to the receiver is
perfectly guaranteed. After getting the reconstructed ${\bf s}^{\prime}$, the
receiver can compute its ${\rm MAC}^{\prime}$ and then be aware of the
integrity of ${\bf s}^{\prime}$ by checking the difference between the
received ${\rm MAC}$ and the computed ${\rm MAC}^{\prime}$.

In light of the central limit theory, the binary sequence $MAC$ follows the
Bernoulli distribution because ${{\bf y}_{MAC}}$ satisfies the sub-Gaussain
distribution. Such a $MAC$ can tolerate some energy-limited noises but some
malicious manipulations.

###  IV-C Access control

In the general outsourced service, the enduser loses the ownership on cloud-
stored data. Hence, the access management is an extremely important work,
which ensure the data confidentiality by restricting access only to the
authorized requesters. After receiving the request for accessing to the
reconstructed data, how to distinguish legitimate access and malicious access
is a significant challenge all the time.

An access passwards generation algorithm based on the random diagonal matrix
${\bf R}_{k_{3}}=diag(r_{1},r_{2},...,r_{N})$ is proposed to restrict
malicious access at a very low complexity. As mentioned previously, the key-
controlled ${\bf R}_{k_{3}}$ plays a pivotal role in the SSR service, whose
diagonal entries are Rademacher random variables. Those diagonal entries cover
the all information which SSR needs such that they are regarded as the shared
secret between the data producer and the cloud provider.

In this paper, those diagonal entries are used to generate the access password
and then it is shared with the legitimate enduser. Let $n$ to be the length of
access password. Here, assume that $n$ is the divisor of $N$. The access
passwords $AP$ can be generated from ${\bf r}$ as follows:

${\bf w}=\frac{1}{2}\left({{{\bf R}_{{k_{3}}}}+{{\bf I}^{N\times
N}}}\right)=diag({w_{1}},{w_{2}},...,{w_{N}})\in\\{0,1\\}$ (15)

where ${\bf I}^{N\times N}$ denotes a $N\times N$ identity matrix. The
acquired ${\bf w}=[{w_{1}},{w_{2}},...,{w_{N}}]$ is a binary sequence. Our
goal is that a $n$-bits random access passward $AP$ is extract from ${\bf w}$
and no information about ${\bf w}$ can be getted from $AP$. Let
$L={N\mathord{\left/{\vphantom{Nn}}\right.\kern-1.2pt}n}$. Two functions $f$
and $g$ are defined as follows:

${\bf p}=f({\bf
w})=f([{w_{1}},{w_{2}},...,{w_{n}}])=[{p_{1}},{p_{2}},...,{p_{n}}]$ (16) ${\bf
q}=f({\bf w})=g([{w_{1}},{w_{2}},...,{w_{n}}])=[{q_{1}},{q_{2}},...,{q_{n}}]$
(17)

where

${p_{i}}=\left\\{{\begin{array}[]{*{20}{c}}{{w_{1}}\otimes{w_{2}}\otimes\cdots\otimes{w_{L}}}\\\
{{p_{i-1}}\otimes{w_{(i-1)L+1}}\otimes{w_{(i-1)L+2}}\otimes\cdots\otimes{w_{iL}}}\end{array}}\right.\begin{array}[]{*{20}{c}}{i=1},\\\
{1<i\leq n},\end{array}$ (18)

${q_{n+1-i}}=\left\\{{\begin{array}[]{*{20}{c}}{{w_{(i-1)L+1}}\otimes{w_{(i-1)L+2}}\otimes\cdots\otimes{w_{nL}}}\\\
{{q_{n-i}}\otimes{w_{(i-1)L+1}}\otimes{w_{(i-1)L+2}}\otimes\cdots\otimes{w_{nL}}}\end{array}\begin{array}[]{*{20}{c}}{i=n},\\\
{1\leq i<n}.\end{array}}\right.$ (19)

the $AP$ is generated by calculating $AP={\bf p}\otimes{\bf q}$, namely

$\begin{array}[]{l}AP={\bf p}\otimes{\bf
q}=[{p_{1}}\otimes{q_{1}},{p_{2}}\otimes{q_{2}},...,{p_{n}}\otimes{q_{n}}]\\\
\end{array}$ (20)

Obviously, the access passwords are $n$ bits and could not leak ${\bf
R}_{k_{3}}$ at least computationally. The access passwords can be generated in
the sensing device and the cloud. The former needs to be secretly sent to the
legitimate end user. After receiving an access request, the cloud manager asks
the requester to submit their access password and then check it. If there
exist any deviation, the requester is viewed as an unauthorized one so this
access is rejected, ortherwise accepted.

##  V The proposed secure outsourcing framework

The section detailedly introduce the proposed secure outsourcing framework
from data encoding and decoding phases. Data encoding is implemented in the
resource-limited sensor but data decoding in both the cloud server and the
terminal device. Note that the fixed orthogonal basis matrix ${\bf\Phi}$ and
hash matrix ${\bf A}_{MAC}$ are kept public. $AP$, $MAC$, and the secret keys
$(k_{1},k_{2},k_{3})$ are secretly transmitted by asymmetric cryptographic
techniques or secure channels.

Figure 2: Privacy-assured outsourcing framework.

###  V-A Data encoding phase

There are three main operations in the encoding phase, including data
acquisition, message authentication codes generation, and access passwords
generation.

_1) Data acquisition:_

  1. Step 1.

Three secret keys $(k_{1},k_{2},k_{3})$ are acquired by key agreement protocol
or a secure channel.

  2. Step 2.

$k_{1}$ and $k_{2}$ are used to generate a random permutation matrix ${{\bf
P}_{k_{1}}}$ and a random diagonal matrix ${{\bf D}_{k_{2}}}$ with non zero
entries, respectively. As a result, the basis matrix can be constructed by
${{\bf\Psi}_{k}}={{\bf\Psi}}{{\bf P}_{k_{1}}}{{\bf D}_{k_{2}}}$.

  3. Step 3.

$k_{3}$ is used to generate a random Rademacher matrix ${\bf R}_{k_{3}}$. As a
result, the sensing matrix ${\bf A}_{k}$ can be constructed by ${\bf
A}_{k}={\bf S}{\bf R}_{k_{3}}{\bf\Psi}$.

  4. Step 4.

The sampling-compression-encryption-hash operation is simultaneously
implemented to capture $\bf y$ and ${\bf y}_{MAC}$ by the following liner
projection:

$\left({\begin{array}[]{*{20}{c}}{\bf y}\\\ {{{\bf
y}_{MAC}}}\end{array}}\right)=\left({\begin{array}[]{*{20}{c}}{{{\bf
A}_{k}}}\\\ {{{\bf
A}_{MAC}}}\end{array}}\right)\cdot{\bf\Psi}_{k}^{-1}\cdot{\bf x}.$ (21)

_2) Message authentication codes generation:_ The $MAC$ of the encrypted
sparse coefficients $\bf s$ are generated by extracting the sign of values as
follows:

$MAC={f_{NRZ}}(sign({{\bf y}_{MAC}})).$ (22)

_3) Access passwords generation:_

  1. Step 1.

The diagonal entries of ${\bf R}_{k_{3}}$ are extracted as a Rademacher random
vector ${\bf r}=[r_{1},r_{2},...,r_{N}]\in\\{-1,1\\}$.

  2. Step 2.

${\bf r}$ is transformed as a symmetric Bernoulli random vector ${\bf
w}=[w_{1},w_{2},...,w_{N}]\in\\{0,1\\}$ by ${\bf w}={\bf r}+1$.

  3. Step 3.

Let $L$ to be step length and
$n=\left\lfloor{{N\mathord{\left/{\vphantom{NL}}\right.\kern-1.2pt}L}}\right\rfloor$.
A $n$-bits binary sequence ${\bf p}=[p_{1},p_{2},...,p_{n}]$ is generated as
follows:

$\begin{array}[]{l}{p_{1}}={w_{1}}\otimes{w_{2}}\otimes\cdots\otimes{w_{L}}\\\
{p_{2}}={p_{1}}\otimes{w_{L+1}}\otimes{w_{L+2}}\otimes\cdots\otimes{w_{2L}}\\\
\vdots\\\
{p_{n}}={p_{n-1}}\otimes{w_{(n-1)L+1}}\otimes{w_{(n-1)L+2}}\otimes\cdots\otimes{w_{nL}}\end{array}$
(23)

  4. Step 4.

A $n$-bits binary sequence ${\bf q}=[q_{1},q_{2},...,q_{n}]$ is generated as
follows:

$\begin{array}[]{l}{q_{1}}={w_{(n-1)L+1}}\otimes{w_{(n-1)L+2}}\otimes\cdots\otimes{w_{nL}}\\\
{q_{2}}={q_{1}}\otimes{w_{(n-2)L+1}}\otimes{w_{(n-2)L+2}}\otimes\cdots\otimes{w_{(n-1)L}}\\\
\vdots\\\
{q_{n}}={q_{n-1}}\otimes{w_{1}}\otimes{w_{2}}\otimes\cdots\otimes{w_{L}}\end{array}$
(24)

  5. Step 5.

A $n$-bits access passwords $AP$ can be generated by $AP={\bf p}\otimes{\bf
q}=[{p_{1}}\otimes{q_{1}},{p_{2}}\otimes{q_{2}},...,{p_{n}}\otimes{q_{n}}]$.
Note that $n$ should not be less than 128 in consideration of security.

###  V-B Data decoding phase

_1) Solving SSR problem:_

  1. Step 1.

The sensing matrix ${\bf A}_{k}={\bf S}{\bf R}_{k_{3}}{\bf\Psi}$ is
constructed through $k$ in the same way as the encoder does.

  2. Step 2.

The privacy-preserving $\bf s$ is got from the measurements $\bf y$ by solving
SSR problem ${\Gamma}_{k}=({\bf A}_{k},{\bf y})$ in the cloud server.

_2) Access control:_

  1. Step 1.

The access passwords $AP$ are generated through $k_{3}$ in the same way as the
encoder does.

  2. Step 2.

After receiving the access request, the correct $AP$ are compared with the
submitted $AP^{\prime}$. If there is no any deviation, this access is
considered to be authorized, otherwise rejected.

_3) Integrity verification:_

  1. Step 1.

After receiving the reconstructed data vector ${\bf s}^{\prime}$, the end-user
computes its $MAC^{\prime}$.

  2. Step 2.

Check $\left\|{MAC-MAC^{\prime}}\right\|\leq\varepsilon$. If this inequality
holds, ${\bf s}^{\prime}$ is viewed to be acceptable, otherwise tampered.

_4) Recovering original data:_

  1. Step 1.

The basis matrix ${{\bf\Psi}_{k}}={{\bf\Psi}}{{\bf P}_{k_{1}}}{{\bf
D}_{k_{2}}}$ is constructed through $k_{1}$ and $k_{2}$ in the same way as the
encoder does.

  2. Step 2.

The original data vector ${\bf x}^{\prime}$ (or called content-preserving
data) is recovered from the verified ${\bf s}^{\prime}$ by computing ${\bf
x}^{\prime}={{\bf\Psi}_{k}}{\bf s}^{\prime}$.

##  VI Experimental results and security analysis

Three kinds of security vulnerabilities are considered in the proposed
outsourcing architecture, including privacy leaking, data tampering, and
malicious access. In this section, we demonstrate the security of the proposed
framework in terms of data confidentiality, data integrity, and access
control. Simulations run in MATLAB R2015B with Core i5-7200U CPU and 4GB RAM.

###  VI-A Confidentiality analysis

Both transmission confidentiality from sensor to cloud and outsourcing
computation privacy are guaranteed by the proposed bi-level encryption
framework, composed of SR-based encryption and CS-based encryption.
Essentially, SR-based encryption is equivalent to performing permutation-
substitution operations in the transformation domain and does not change the
sparsity. For more intuitive experimental results, its encrypted results are
exhibited in the time domain. CS-based encryption aims at data obfuscation
while randomly subsampling. Such a product cipher is perfectly built in the
sensing player at almost zero cost and provides the captured data with innate
confidentiality. Here, we investigate its privacy-preserving capability by
simulation experiments and theoretical analysis below.

Compressive ratio $CR$ is set to 0.5. The four plain images are encrypted
through SR, then encrypted through CS, and finally decrypted as shown in Fig.
3. The finally encrypted image containing less data are transmitted from the
sensor and the cloud. The halfway encrypted image contains all available
information for the cloud. Obviously, both the halfway encrypted images and
the finally encrypted images do not leak any visually meaningful information.
Note that the halfway encrypted image leaks the sparsity of the original
signal. The visual secrecy provided by SR-based encryption is acceptable in
the face of the semi-trusted cloud but not enough to guarantee secure
transmission over the public channels. Such a vulnerability occurred in [18,
19, 20]. Next, we consider some common attack scenarios in the process of
transmission, including brute-force attack, ciphertext-only attack (COA), and
plaintext attack.

_1) Brute-force attack:_ As mentioned previously, both the sensing matrix and
the basis matrix are controlled by several keys, namely
${{\bf\Psi}_{k}}={{\bf\Psi}}{{\bf P}_{k_{1}}}{{\bf D}_{k_{2}}}$, ${\bf
A}_{k}={\bf S}{\bf R}_{k_{3}}{\bf\Psi}$. The whole encoding process depends on
three random matrix $({{\bf P}_{k_{1}}},{{\bf D}_{k_{2}}},{\bf R}_{k_{3}})$.
Hence, the potential case of brute-force attack is that the attacker tries to
get the real matrices $({{\bf P}_{k_{1}}},{{\bf D}_{k_{2}}},{\bf R}_{k_{3}})$
by guessing or searching the whole key space of $(k_{1},k_{2},k_{3})$
exhaustively and then uses them for decoding.

Let $Pro(\cdot)$ denotes the probability of a successful trial. It is not hard
to calculate $Pro({{\bf
P}_{k_{1}}})={1\mathord{\left/{\vphantom{1N}}\right.\kern-1.2pt}N}!$,
$Pro({{\bf
D}_{k_{2}}})={1\mathord{\left/{\vphantom{1{{2^{N}}}}}\right.\kern-1.2pt}{{f^{N}}}}$,
$Pro({{\bf
R}_{k_{3}}})={1\mathord{\left/{\vphantom{1{{2^{N}}}}}\right.\kern-1.2pt}{{2^{N}}}}$.
So the probability of successful guessing the three matrices is $Pro({{\bf
P}_{k_{1}}},{{\bf D}_{k_{2}}},{\bf
R}_{k_{3}})={1\mathord{\left/{\vphantom{1{(N!\cdot{f^{N}}\cdot{2^{N}}}}}\right.\kern-1.2pt}{(N!\cdot{f^{N}}\cdot{2^{N}}}})$.
For a $N$-conditionality multimedia data, $N$ is often enough large to make
this brute-force attack infeasible, namely $Pro({{\bf P}_{k_{1}}},{{\bf
D}_{k_{2}}},{\bf R}_{k_{3}})\\\ \approx 0$.

_2) Cipher text-only attack:_ A recent work pointed out that CS measurements
acquired by a structurally subsampled sensing matrix are said to be
indistinguishable as long as the original signal has constant energy[32]. Such
a CS-based cryptosystem can achieve the asymptotic spherical secrecy[33]. It
means that $\bf y$ only leak the energy information of $\bf x$ when a
structurally subsampled sensing matrix is employed, rooted in the energy-
preserving constraint of RIP condition. However, a non-normalized random
matrix ${{\bf D}_{k_{2}}}$ is imported to CS-based encoding in the proposed
framework. It can break the energy conservation between $\bf x$ and $\bf y$
but does not affect the exact reconstruction such that a higher secrecy than
asymptotic spherical secrecy is achieved in the proposed framework.

###### Theorem 5

[32] A CS-based cryptosystem ${\bf y}={\bf A}{\bf x}$ can be said to achieve
asymptotic spherical secrecy as long a $\bf A$ is a structurally subsampled
sensing matrix, where $\bf x$ and $\bf y$ denote the spare plaintext and the
corresponding ciphertext, respectively.

###### Deifinition 2

(Asymptotic spherical secrecy[33]) Define ${\bf
x}=[{x_{1}},{x_{2}},{x_{3}},...,{x_{n}}]$ as a plaintext sequence and ${\bf
y}$ as the corresponding ciphertext sequence. Assume that the energy of the
plaintext ${\varepsilon_{\bf x}}$ is finite. Such a cryptosystem is said to
achieve asymptotic spherical secrecy if ${f_{{\bf y}|{\bf
x}}}(x,y)\mathop{\to}\limits_{\mathcal{D}}{f_{{\bf y}|{\varepsilon_{\bf
x}}}}(y)$, where ${\mathcal{D}}$ denotes the convergence in the distribution
as $n\to\infty$.

_2) Plaintext attack:_ From the perspective of stealing the complete
plaintext, It is widely accepted that CS-based basic cryptosystem is immune to
brute-force attack and COA. Special attention is paid to the plaintext
attacks. Plaintext attack means that abundant plaintext-ciphertext pairs have
been captured. When the measurement matrix $\bf\Phi$ is fixed, the attacker
can figure out $bf\Phi$ from $N$ independent plaintext-ciphertext pairs and
further illegally recover the current plaintext by leveraging $\bf\Phi$. To
resist plaintext attack, basic CS cipher is forced to renew $\bf\Phi$
frequently, even in a one-time sensing (OTS) model.

In the proposed framework, the measurement matrix is $\bf\Phi={\bf S}{\bf
R}_{k_{3}}{\bf\Psi}({{\bf\Psi}}{{\bf P}_{k_{1}}}{{\bf D}_{k_{2}}})^{-1}$.
${\bf S}{\bf R}_{k_{3}}{\bf\Psi}$ has been demonstrated to be a RIP matrix.
${{\bf D}_{k_{2}}}$ breaks the energy-preserving guarantee of RIP such that
$\bf\Phi$ is a non-RIP matrix. In a multi-time sensing (MTS) model, $\bf\Phi$
can also be calculated from $N$ independent plaintext-ciphertext pairs. The
difficulty of plaintext attack depends on the complexity of getting $\bf D$ by
decomposing $\bf\Phi$ or guessing. No matter which approaches, it seems to be
an impossible task in polynomial time.

Therefore, without updating the key $(k_{1},k_{2},k_{3})$ frequently, the
proposed framework can provide outsourced data with a high level of
confidentiality against brute-force attack, COA, and plaintext attack.

(a) Plain images (b) Halfway encrypted images (c) Finally encrypted images (d)
Decrypted images

Figure 3: Encrypted and decrypted results for four plain images in $CR=0.5$.

###  VI-B Integrity analysis

A noise-robust but tamper-resistant integrity verification approach is
involved in the proposed framework. Considering that both transmission and
outsourcing computations carry out in the ciphertext domain, the target of
tampering the plaintext data $\bf x$ seems to be achieved only in two cases.
One that polluting measurements $\bf y$ makes the reconstruction process fail
or the reconstructed data beyond recognition. The other is that $\bf y$ is
replaced with anther fake containing entirely different information, based on
the worst scenario where the eavesdropper has grasped the practical
measurement matrix $\bf\Phi$. Therefore, we discuss four situations as
described in Fig. 4, namely no noise, normal noise, malicious pollution, and
fake image. It is no surprise that $MAC$ generated from the recovered data are
different from one generated from the original data. The bit error rate (BER)
of $MAC$ is used to estimate the distortion degree of $MAC$. The experimental
results are shown in Fig. 5. There is almost no change in $MAC$ under no noise
circumstance. The normal noise makes a slight difference to $MAC$ but
malicious noise dose the opposite. The fake data can also change $MAC$ to some
extent. Obviously, data pollution and replacement behaviors can be easily
detected when the threshold value of BER is set in the interval $(0.1,0.15)$.
The performance of such a data integrity verification approach is affected by
CR, noise, reconstruction algorithm, etc. Note that some small-region data
tamper may be viewed as the normal noise.

Figure 4: The reconstructed results in four different situations. (a) no
noise; (b) normal noise; (c) malicious pollution; (c) fake image. Figure 5:
BER versus CR in different situations.

###  VI-C Access control analysis

To restrict the malicious access, a lightweight access passwords $AP$
generation algorithm is proposed in Section 4, which is based on the prior
knowledge about sensing matrix. The end-user launches an access request by
submitting $AP$ and the cloud manager verifies the validity of this access by
checking $AP$. The transmission security of $AP$ depends on the key agreement
protocols, which is outside the scope of this paper. Taking security into
consideration, we suggest that $AP$ should be at least 128-bits. Here, we test
the sensitivity of access passwords generation algorithm. Essentially, the
proposed algorithm is to transform a length-variable Bernoulli random sequence
to an unpredictable 128-bits binary sequence. Three kinds of data
manipulations are performed on the input sequence, including inverting some
elements, randomly permuting two elements, slightly key perturbation. As shown
in Table 1, the output AP would be radically changed no matter which
manipulation is performed. It means that each encoding would generate an
exclusive access password with high probability. Only by submitting a correct
access password can the access requester be considered as an authorized user.

TABLE I: Access password test Manipulation types |  |  Inverse first element |  Inverse middle element |  Inverse last element |  Permutate two elements |  Key perturbation
---|---|---|---|---|---|---
Access passwords |  6BF27A9B B956A7BA 5FDF81DB C394F691 |  21DFF11E CF2EC346 85EE6CA9 A841AA7E |  38D0BFD1 109A1BF2 9D09D896 CE0340C6 |  D3D1AF9A 6BE27A9B B956B7BA 7F4E83CB |  E2EA6787 89AAFA71 06D72E17 A12C67DC |  725E7524 3C61D8ED F46B4103 09036C4A
Change rate |  |  95.14% |  97.73% |  95.62% |  94.35% |  96.81%

##  VII Conclusion

This paper considers such a scenario in which both sensing device and terminal
device are resource-constrained but the cloud is resource-abundant, which
generally appears in the current IoT applications. With the boom of multimedia
big data, it is an predictable trend to transform computational complexity
form the local side to the cloud side. In this paper, a privacy-assured
outsourcing scheme is proposed for multimedia, in which CS plays a pivotal
role. CS’s low-complexity encoding is carried forward in the resource-
constrained sensing device and high-complexity decoding is left to the cloud.
More importantly, several security issues are effectively solved in CS
framework.

A CS-based product cipher is proposed to provide the outsourced data with
secrecy guarantee in the process of transmission and computation. Security
analyses indicate that the captured data are able to resist brute-force
attack, COA, and plaintext attack. Considering the malicious access attack to
cloud, an access password generation algorithm is proposed, which is based on
the prior knowledge of the structurally subsampled sensing matrix. In
addition, data integrity guarantee is also acquired by a CS-based message
authentication codes generator. Differing from most sensitive hash algorithm,
CS-based message authentication codes can not only resist noise and tamper
attacks but tolerate energy-finite noise to some extent.

## References

  * [1] M. Lustig, D. L. Donoho, J. M. Santos, J. M. Pauly, Compressed sensing mri, IEEE signal proc. mag. 25 (2) (2008) 72–82.
  * [2] C. R. Berger, Z. Wang, J. Huang, S. Zhou, Application of compressive sensing to sparse channel estimation, IEEE Commun. Maga. 48 (11) (2010) 164–174.
  * [3] S. Li, L. Da Xu, X. Wang, Compressed sensing signal and data acquisition in wireless sensor networks and internet of things, IEEE Trans. Ind. Inform. 9 (4) (2012) 2177–2186.
  * [4] M. J. Atallah, K. N. Pantazopoulos, J. R. Rice, E. E. Spafford, Secure outsourcing of scientific computations, in: adv. comput., Vol. 54, Elsevier, 2002, pp. 215–272.
  * [5] S. Hohenberger, A. Lysyanskaya, How to securely outsource cryptographic computations, in: Proc. Theory Cryptogr. Conf. (TCC), Springer, 2005, pp. 264–282.
  * [6] M. J. Atallah, K. B. Frikken, Securely outsourcing linear algebra computations, in: Proc. 5th ACM Symp. Inf. Comput. Commun. Security (ASIACCS), 2010, pp. 48–59.
  * [7] F. Kerschbaum, Outsourced private set intersection using homomorphic encryption, in: Proc. ACM Symp. Inf. Comput. Commun. Security (ASIACCS), 2012, pp. 85–86.
  * [8] M. Ibtihal, D. El Ouadghiri, N. Hassan, Homomorphic encryption as a service for outsourced images in mobile cloud computing environment, Int. J. Cloud Appl. Comput. 7 (2017) 27–40.
  * [9] P. Li, J. Li, Z. Huang, C.-Z. Gao, W.-B. Chen, K. Chen, Privacy-preserving outsourced classification in cloud computing, Cluster Comput. 21 (1) (2018) 277–286.
  * [10] X. Lei, X. Liao, T. Huang, H. Li, C. Hu, Outsourcing large matrix inversion computation to a public cloud, IEEE Trans. cloud comput. 1 (1) (2013) 1–1.
  * [11] X. Lei, X. Liao, T. Huang, F. Heriniaina, Achieving security, robust cheating resistance, and high-efficiency for outsourcing large matrix multiplication computation to a malicious cloud, Inf. sci. 280 (2014) 205–217.
  * [12] X. Lei, X. Liao, T. Huang, H. Li, Cloud computing service: The caseof large matrix determinant computation, IEEE Trans. Serv. Comput. 8 (5) (2014) 688–700.
  * [13] C. Wang, K. Ren, J. Wang, K. M. R. Urs, Harnessing the cloud for securely solving large-scale systems of linear equations, in: 31st Int. Conf. Distrib. Comput. Syst., 2011, pp. 549–558.
  * [14] C. Wang, K. Ren, J. Wang, Q. Wang, Harnessing the cloud for securely outsourcing large-scale systems of linear equations, IEEE Trans. Parallel Distrib. Syst. 24 (6) (2013) 1172–1181.
  * [15] C. Wang, K. Ren, J. Wang, Secure and practical outsourcing of linear programming in cloud computing, in: Proc. IEEE Computer Commun, 2011, pp. 820–828.
  * [16] C. Wang, K. Ren, J. Wang, Secure optimization computation outsourcing in cloud computing: A case study of linear programming, IEEE Transactions on Computers 65 (1) (2016) 216–229.
  * [17] C. Wang, B. Zhang, K. Ren, J. M. Roveda, Privacy-assured outsourcing of image reconstruction service in cloud, IEEE Trans. Emerg. Top. Comput. 1 (1) (2013) 166–177.
  * [18] Y. Zhang, J. Zhou, L. Y. Zhang, F. Chen, X. Lei, Support-set-assured parallel outsourcing of sparse reconstruction service for compressive sensing in multi-clouds, in: Int. Symp. Secur. Priv. Soc. Netw. Big Data (SocialSec), IEEE, 2015, pp. 1–6.
  * [19] G. Hu, D. Xiao, T. Xiang, S. Bai, Y. Zhang, A compressive sensing based privacy preserving outsourcing of image storage and identity authentication service in cloud, Inf. Sci. 387 (2017) 132–145.
  * [20] Y. Zhang, Y. Xiang, L. Y. Zhang, L.-X. Yang, J. Zhou, Efficiently and securely outsourcing compressed sensing reconstruction to a cloud, Inf. Sci. 496 (2019) 150–160.
  * [21] L. Y. Zhang, K. W. Wong, Y. Zhang, J. Zhou, Bi-level protected compressive sampling, IEEE Trans. Multimedia 18 (9) (2016) 1720–1732.
  * [22] Y. Zhang, B. Xu, N. Zhou, A novel image compression–encryption hybrid algorithm based on the analysis sparse representation, Opt. Commun. 392 (2017) 223–233.
  * [23] N. Zhou, H. Jiang, L. Gong, X. Xie, Double-image compression and encryption algorithm based on co-sparse representation and random pixel exchanging, Opt. Lasers Eng. 110 (2018) 72–79.
  * [24] D. L. Donoho, Compressed sensing, IEEE Trans. Inf. Theory 52 (4) (2006) 1289–1306.
  * [25] E. J. Candes, M. B. Wakin, An introduction to compressive sampling, IEEE Signal Process. Mag. 25 (2) (2008) 21–30.
  * [26] E. J. Candes, J. Romberg, T. Tao, Robust uncertainty principles: exact signal reconstruction from highly incomplete frequency information, IEEE Trans. Inf. Theory 52 (2) (2006) 489–509.
  * [27] E. J. Candes, T. Tao, Near-optimal signal recovery from random projections: universal encoding strategies?, IEEE Trans. Infor. Theory 52 (12) (2006) 5406–5425.
  * [28] R. Baraniuk, M. Davenport, R. DeVore, M. Wakin, A simple proof of the restricted isometry property for random matrices, Constr. Approx. 28 (3) (2008) 253–263.
  * [29] R. Kueng, D. Gross, Ripless compressed sensing from anisotropic measurements, Linear Algebra App. 441 (2014) 110–123.
  * [30] W. U. Bajwa, A. M. Sayeed, R. Nowak, A restricted isometry property for structurally-subsampled unitary matrices, in: 47th annu Allerton conf commun control, and comput (Allerton), IEEE, 2009, pp. 1005–1012.
  * [31] M. F. Duarte, Y. C. Eldar, Structured compressed sensing: From theory to applications, IEEE Trans. signal proces. 59 (9) (2011) 4053–4085.
  * [32] W. Cho, N. Y. Yu, Secure and efficient compressed sensing-based encryption with sparse matrices, IEEE Trans. Inf. Foren. Sec. 15 (2019) 1999–2011.
  * [33] V. Cambareri, M. Mangia, F. Pareschi, R. Rovatti, G. Setti, Low-complexity multiclass encryption by compressed sensing, IEEE Trans. Signal Process. 63 (9) (2015) 2183–2195.


# Explosive behaviour in networks of Winfree oscillators

Shawn Means  School of Natural and Mathematical Sciences, Massey University,
Private Bag 102-904 North Shore Mail Centre, Auckland 0745, New Zealand
<EMAIL_ADDRESS>and  Carlo R. Laing  School of Natural and Mathematical
Sciences, Massey University, Private Bag 102-904 North Shore Mail Centre,
Auckland 0745, New Zealand<EMAIL_ADDRESS>
###### Abstract.

We consider directed networks of Winfree oscillators with power law
distributed in- and out-degree distributions. Gaussian and power law
distributed intrinsic frequencies are considered, and these frequencies are
positively correlated with oscillators’ in-degrees. The Ott/Antonsen ansatz is
used to derive degree-based mean field equations for the expected dynamics of
networks, and these are numerically analysed. In a variety of cases
“explosive” transitions between either two different steady states or between
a steady state and a periodic solution are found, and these transitions are
explained using bifurcation theory.

##  1\. Introduction

The study of networks is driven by the many real-world systems which can be
represented in this way [23], such as networks of neuronal [28], cardiac [4]
or smooth muscle [40] cells. In networks of oscillators it is often the
conditions under which the oscillators synchronise which are of interest [29,
36]. For networks of nonidentical oscillators there is often a continuous
transition from asynchrony to higher and higher levels of synchrony as the
strength of coupling between oscillators is increased [12, 36]. But recently a
different form of transition known as “explosive synchrony” (ES) has been
observed. Here, there is a discontinuous increase in network coherence as the
coupling strength between oscillators is increased (see [6] for a review).
This behaviour was first reported by [8] in an undirected network of Kuramoto
oscillators [12] for which the intrinsic frequency of an oscillator was equal
to its degree, the degrees having an inverse power law distribution. The
authors emphasised that the explosive transition was due to there being a
relationship between a local network property (an oscillator’s degree) and a
property of an oscillator (its intrinsic frequency) — in this case a positive
correlation between the two. Since the publication of [8] there have been many
studies of ES [17, 2, 42, 41].

As recently noted by [11], ES is not actually an unexpected phenomenon, since
if there is a continuous transition to synchrony through either a
transcritical or pitchfork bifurcation as, say, the strength of coupling
between oscillators is increased, then by varying a second parameter
sufficiently the criticality of the bifurcation will change. This change leads
to bistability and a region of hysteresis as the original parameter is varied
(i.e. an “explosive” transition).

Many studies of ES use Kuramoto phase oscillators which are coupled through
phase differences. However, more realistic oscillators are not coupled in this
way, instead having interactions which depend on the states of the two
oscillators which are coupled [37, 21, 3, 19]. The oscillator we study – the
Winfree model [39] – is a phase oscillator similar in form to the Kuramoto
oscillator but is not described with phase differences, instead involving a
“phase response curve” and a pulsatile function of phase. We are not aware of
any reports of ES in networks of Winfree oscillators, and as one of our
contributions we report here our observations of this phenomenon in our
numerical explorations of networks of Winfree oscillators in a variety of
network configurations.

Most previous studies of networks of Winfree oscillators have focused on the
all-to-all coupled case [1, 26, 10, 7, 27], with an exception being [16]. Many
previous studies of ES have considered undirected networks, where oscillators
are either coupled or not. However, we consider directed networks – as seen in
neuronal systems – where one oscillator influences other oscillators coupled
downstream but without necessarily a reciprocal influence back upstream.

The bistability that often accompanies ES is of interest as it allows a
network to have more than one stable state for a given set of parameters, with
the current state often determined by the network’s history. A prominent
example of such bistability is in networks which exhibit Up and Down states
[38, 30], with transitions between these states driven by either fluctuations
or noise of some form, or by transient inputs. The abrupt transitions between
states may also provide a model for the conscious-unconscious transition when
awaking from anesthesia [35].

The structure of the paper is as follows. In Sec. 2 we present the network we
consider and the degree-based mean field description of its dynamics, derived
using the Ott/Antonsen ansatz. Specific cases of independent in- and out-
degrees, and equal in- and out-degrees, are considered. In Sec. 3 we consider
a Gaussian distribution of intrinsic frequencies which are positively
correlated with oscillators’ in-degrees. In Sec. 4 we consider an inverse
power law distribution of intrinsic frequencies and again correlate them with
oscillators’ in-degrees. We consider the cases of independent in- and out-
degrees, and equal in- and out-degrees. A variety of “explosive” transitions
are observed and explained using simple bifurcation theory. We conclude in
Sec. 5.

##  2\. Model

The model describing a directed network of Winfree oscillators is as presented
in [16, 26]:

$\frac{d\theta_{j}}{dt}=\omega_{j}+U(\theta_{j})\frac{\epsilon}{\langle
k\rangle}\sum_{n=1}^{N}A_{jn}T(\theta_{n});\qquad j=1,2\dots N$ (1)

where $\epsilon$ is the strength of coupling, $\langle k\rangle$ is the mean
degree of the network and the connectivity is given by the adjacency matrix
$A$ with $A_{jn}=1$ if oscillator $n$ connects to oscillator $j$ and zero
otherwise. The $\omega_{j}$ are chosen from a Lorentzian distribution with
centre and half-width at half-maximum to be discussed below. The phase
response curve [32, 22] is chosen to be $U(\theta)=-\sin{\theta}$ and the
pulsatile function $T$ which has a maximum at $\theta=0$ is

$T(\theta)=\frac{8}{35}(1+\cos{\theta})^{4}$ (2)

The in-degree of oscillator $j$ is

$k_{in,j}=\sum_{n=1}^{N}A_{jn}$ (3)

and the out-degree of oscillator $n$ is

$k_{out,n}=\sum_{j=1}^{N}A_{jn}$ (4)

We consider large networks where all in- and out-degrees are large and assume
that the network can be characterised by its degree distribution $P({\bf k})$
where ${\bf k}=(k_{in},k_{out})$, and an assortativity function $a({\bf
k}^{\prime}\to{\bf k})$ giving the probability that an oscillator with degree
${\bf k}^{\prime}$ connects to one with degree ${\bf k}$, given that such
oscillators exist. $P({\bf k})$ is normalised such that $\sum_{{\bf k}}P({\bf
k})=N$ and we choose the marginal distributions of the in-degrees and the out-
degrees to be equal.

Using the theory in [16] (or see [15, 5] for similar derivations) one can show
that the long-time dynamics of the network is described by

$\frac{\partial b({\bf k},t)}{\partial t}=\frac{\epsilon R({\bf
k},t)}{2}+[i\omega_{0}({\bf k})-\Delta({\bf k})]b({\bf k},t)-\frac{\epsilon
R({\bf k},t)}{2}[b({\bf k},t)]^{2}$ (5)

where $\omega_{0}({\bf k})$ and $\Delta({\bf k})$ are the centre and half-
width at half-maximum, respectively, of the Lorentzian distribution from which
the values of $\omega_{j}$ for oscillators with degree ${\bf k}$ are chosen:

$g(\omega({\bf k}))=\frac{\Delta({\bf k})/\pi}{[\omega({\bf
k})-\omega_{0}({\bf k})]^{2}+\Delta^{2}({\bf k})}$ (6)

The variable

$b({\bf k},t)=\int_{-\infty}^{\infty}\int_{0}^{2\pi}f(\theta,\omega|{\bf
k},t)e^{-i\theta}d\theta\ d\omega$ (7)

is the complex-valued order parameter for oscillators with degree ${\bf k}$,
where $f(\theta,\omega|{\bf k},t)d\theta d\omega$ is the probability that an
oscillator with degree ${\bf k}$ has phase in $[\theta,\theta+d\theta]$ and
frequency in $[\omega,\omega+d\omega]$ at time $t$. $R({\bf k},t)$ is given by

$R({\bf k},t)=\frac{1}{\langle k\rangle}\sum_{{\bf k}^{\prime}}P({\bf
k}^{\prime})a({\bf k}^{\prime}\to{\bf k})G({\bf k}^{\prime},t)$ (8)

where

$\displaystyle G({\bf k},t)$ $\displaystyle=$ $\displaystyle
1+\frac{4\left[b({\bf k},t)+\overline{b}({\bf
k},t)\right]}{5}+\frac{2\left[b^{2}({\bf k},t)+\overline{b}^{2}({\bf
k},t)\right]}{5}$ (9) $\displaystyle+\frac{4\left[b^{3}({\bf
k},t)+\overline{b}^{3}({\bf k},t)\right]}{35}+\frac{b^{4}({\bf
k},t)+\overline{b}^{4}({\bf k},t)}{70}$

where overline indicates the complex conjugate. The form of $G$ is determined
by the function $T(\theta)$. This derivation of a degree-dependent mean field
description uses the Ott/Antonsen ansatz [24, 25]. A crucial ingredient for
the use of this ansatz is the sinusoidal form of the function $U(\theta)$.

We are interested in the case of neutral assortativity, for which [31]

$a({\bf k}^{\prime}\to{\bf k})=\frac{k_{out}^{\prime}k_{in}}{N\langle
k\rangle}$ (10)

and cases where either the in- and out-degree of an oscillator are
independent, or they are equal (and thus perfectly correlated). Writing
$P(k_{in}^{\prime},k_{out}^{\prime})$ instead of $P({\bf k}^{\prime})/N$ we
have

$R(k_{in},k_{out})=\frac{k_{in}}{\langle
k\rangle^{2}}\sum_{k_{in}^{\prime}}\sum_{k_{out}^{\prime}}P(k_{in}^{\prime},k_{out}^{\prime})k_{out}^{\prime}G(k_{in}^{\prime},k_{out}^{\prime},t)$
(11)

which is clearly independent of $k_{out}$. Thus both $b$ and $G$ must also be
independent of $k_{out}$ and we have

$R(k_{in},t)=\frac{k_{in}}{\langle
k\rangle^{2}}\sum_{k_{in}^{\prime}}Q(k_{in}^{\prime})G(k_{in}^{\prime},t)$
(12)

where

$Q(k_{in}^{\prime})=\sum_{k_{out}^{\prime}}P(k_{in}^{\prime},k_{out}^{\prime})k_{out}^{\prime}.$
(13)

For the case of the in- and out-degrees of an oscillator being independent,
$P(k_{in}^{\prime},k_{out}^{\prime})$ factorises:
$P(k_{in}^{\prime},k_{out}^{\prime})=p(k_{in})p(k_{out}^{\prime})$ where $p$
is the marginal distribution of either the in- or out-degrees. In this case,

$R(k_{in},t)=\frac{k_{in}}{\langle
k\rangle}\sum_{k_{in}^{\prime}}p(k_{in}^{\prime})G(k_{in}^{\prime},t).$ (14)

Alternatively, if $k_{in}=k_{out}$ for each oscillator (i.e. each oscillator
has the same in- and out-degree) then
$Q(k_{in}^{\prime})=k_{in}^{\prime}p(k_{in}^{\prime})$ and

$R(k_{in},t)=\frac{k_{in}}{\langle
k\rangle^{2}}\sum_{k_{in}^{\prime}}p(k_{in}^{\prime})k_{in}^{\prime}G(k_{in}^{\prime},t).$
(15)

In either case the “input” to an oscillator with in-degree $k_{in}$ is
proportional to its in-degree, i.e. $k_{in}$.

One could use this theory to investigate the effects of choosing different
degree distributions $p$ as in [14] but instead we fix $p$ as a power law
distribution and look at the effects of correlating an oscillator’s frequency
with its degree. Specifically, we make $\Delta$ independent of degree but
choose $\omega_{0}$ to be a function of $k_{in}$ only. (One could make
$\omega_{0}$ a function of $k_{out}$ as well, or instead, but doing so
increases the computational cost, and the phenomena we are interested in occur
when $\omega_{0}$ is a function of $k_{in}$ only.) Thus the equations of
interest are

$\frac{\partial b(k_{in},t)}{\partial t}=\frac{\epsilon
R(k_{in},t)}{2}+[i\omega_{0}(k_{in})-\Delta]b(k_{in},t)-\frac{\epsilon
R(k_{in},t)}{2}[b(k_{in},t)]^{2}$ (16)

where $m\leq k_{in}\leq M$, where $m$ and $M$ are the minimum and maximum in-
and out-degrees respectively, $R$ is given by either (14) or (15) and $G$ is
given by (9) but with only dependence on in-degree.

Note that in previous work [16] we used similar equations as those above to
investigate the effects of varying the correlation between in- and out-degrees
of Winfree oscillators, and of parameter- and degree-based assortativity. We
also examined the effects of correlating oscillators’ intrinsic frequencies
with their in-degrees, but did not consider varying the coupling strength.
Also, there we considered uniform degree distributions rather than power law,
as we do here and as many other researchers do.

Below we will use the mean firing frequency of a network to describe its
behaviour, so now derive an expression for that. The dynamics of an oscillator
with in-degree $k_{in}$ is

$\frac{d\theta}{dt}=\omega(k_{in})-\epsilon R(k_{in})\sin{\theta}.$ (17)

If $|\omega(k_{in})|>\epsilon R(k_{in})$ then the oscillator will fire
periodically with frequency

$\frac{\sqrt{\omega^{2}(k_{in})-\epsilon^{2}R^{2}(k_{in})}}{2\pi}.$ (18)

Thus the expected frequency for oscillators with in-degree $k_{in}$ is

$\displaystyle f(k_{in})$
$\displaystyle=\frac{1}{2\pi}\mbox{Re}\left(\int_{-\infty}^{\infty}g(\omega(k_{in}))\sqrt{\omega^{2}(k_{in})-\epsilon^{2}R^{2}(k_{in})}\
d\omega(k_{in})\right)$
$\displaystyle=\frac{1}{2\pi}\sqrt{\frac{\omega_{0}^{2}(k_{in})-\Delta^{2}-\epsilon^{2}R^{2}(k_{in})+\sqrt{[\omega_{0}^{2}(k_{in})-\Delta^{2}-\epsilon^{2}R^{2}(k_{in})]^{2}+4\Delta^{2}\omega_{0}^{2}(k_{in})}}{2}}$
(19)

where $g$ is the Lorentzian (6) and thus the mean firing rate over the network
is

$F=\frac{1}{2\pi}\sum_{k_{in}}p(k_{in})f(k_{in})$ (20)

For comparison with our results below we briefly describe the dynamics of a
fully connected network of Winfree oscillators. Two types of behaviour are
typically seen in such a network: synchronous and asynchronous [26, 7],
although the fraction of oscillators actually oscillating varies in different
asynchronous states. Increasing $\epsilon$ tends to destroy synchronous
behaviour through a saddle-node-on-invariant-circle (SNIC) bifurcation, as
many of the oscillators “lock” to an approximate fixed point. For moderate
$\epsilon$ increasing the spread of intrinsic frequencies tends to destroy
synchronous behaviour through a supercritical Hopf bifurcation, as the
oscillators become too dissimilar in frequency to synchronise [26]. Below we
will see a wider variety of bifurcations resulting from the networks’
structure.

##  3\. Gaussian frequency distribution

We choose the degree distribution $p(k)$ to be a truncated power law
distribution with exponent $-3$, as many others have done when studying ES
[18, 8]:

$p(k)=\begin{cases}a/k^{3}&m\leq k\leq M\\\ 0&\mbox{ otherwise}\end{cases}$
(21)

where $a$ is a normalisation such that

$\sum_{k=m}^{M}\frac{a}{k^{3}}=1.$ (22)

Since the degrees are all large (i.e. $1\ll m$) we treat $k$ as a continuous
variable and approximate the sum in (22) by an integral, giving
$a=2m^{2}M^{2}/(M^{2}-m^{2})$. The cumulative distribution function for $k$ is

$\widehat{p}(k)=\int_{m}^{k}p(s)\
ds=\frac{a}{2}\left(\frac{1}{m^{2}}-\frac{1}{k^{2}}\right).$ (23)

We need to specify the dependence of $\omega_{0}$ on $k_{in}$. In this section
we consider the case where for a particular realisation of the network, we
randomly choose the target frequencies from a Gaussian distribution with mean
1 and standard deviation $\sigma$, then assign the smallest target frequency
to the oscillator with smallest in-degree, the second smallest target
frequency to the oscillator with second smallest in-degree, etc. (For
oscillators with equal in-degree, they are ranked in random order.) For
oscillator $j$ the actual $\omega_{j}$ is then chosen from a Lorentzian with
centre equal to the target frequency for oscillator $j$ and half-width at
half-maximum (HWHM) equal to $\Delta$, similar to the scheme in [34].

In this case we have

$\omega_{0}(k)=\widehat{q}^{-1}(\widehat{p}(k))$ (24)

where $\widehat{q}$ is the cumulative distribution function of the appropriate
Gaussian distribution, i.e.

$\widehat{q}(\omega)=\frac{1}{2}\left[1+\mbox{erf}\left(\frac{\omega-1}{\sqrt{2}\sigma}\right)\right]$
(25)

A demonstration of this is shown in Fig. 1 where we create a directed network
with $N=2000$ nodes using the configuration model [23], then randomly chose
$N$ target frequencies from a Gaussian distribution with mean 1 and standard
deviation $0.01$, then associated the smallest frequency with the node with
smallest in-degree etc. These values are shown as dots in panel (a), along
with the theoretical relationship (24). Panel (b) shows the actual values of
the $\omega_{j}$ (dots), taken from a Lorentzian with centre equal to the
values in panel (a) and HWHM 0.00005.

Figure 1. (a): solid line: (24); dots: target frequencies and degrees for one
realisation of a network. (b): solid line: (24); dots: actual frequencies and
degrees for one realisation of a network. Parameters:
$m=100,M=400,N=2000,\sigma=0.01,\Delta=0.00005$.

Clearly there is a positive correlation between $k_{in,j}$ and $\omega_{j}$,
but the Pearson correlation coefficient between them, defined by

$\rho_{k\omega}=\frac{\sum_{j=1}^{N}(k_{in,j}-\langle
k\rangle)(\omega_{j}-\bar{\omega})}{\sqrt{\sum_{j=1}^{N}(k_{in,j}-\langle
k\rangle)^{2}}\sqrt{\sum_{j=1}^{N}(\omega_{j}-\bar{\omega})^{2}}}$ (26)

where $\bar{\omega}$ is the mean of the $\omega_{j}$, is less than 1 and
(using this construction) cannot be systematically varied, as was possible in
[16, 15]. Note that the idea of not having a perfect match between an
oscillator’s degree and a prescribed frequency (as we have here for
$\Delta\neq 0$) was discussed in [33], where it was shown that such a mismatch
actually created ES.

To investigate the influence of varying the correlation between $k_{in,j}$ and
$\omega_{j}$, we created a sequence of the appropriate degrees, chose values
of $\omega$ from a Gaussian distribution and randomly paired them. We then
repeatedly performed Monte Carlo swaps of $\omega$ values; potential swaps
were accepted if they increased the Pearson correlation coefficient toward a
target value – typically a positive number. This approach required substantial
effort to satisfy high correlations $(\geq 0.9)$, since a sequence of $\omega$
randomly paired with degree sequences typically exhibited no correlation.
Alternatively, we maximised the correlation between the $k_{in,j}$ and
$\omega_{j}$ by initially sorting both sequences as above. Aligning maximal-
minimal values then produced the highest possible correlation coefficient for
given sequences which we then reduced using Monte Carlo swaps, with swaps
accepted if they pushed the correlation value toward a target value. Sequences
of degrees and frequencies generated with this approach were assembled into
adjacency matrices utilising our network assembly scheme, the “permutation
method”, presented previously [20]. Note that with this approach of
constructing networks the Ott/Antonsen approach cannot be used, and we must
simulate the resulting networks to determine their behaviour.

###  3.1. Results

We numerically investigate (16) with (14). We evaluate functions at all
integer in-degrees satisfying $m+1\leq k_{in}\leq M-1$ to avoid the
singularities in $\widehat{q}^{-1}$ when its argument is either 0 or 1,
resulting in a moderately large set of ordinary differential equations. We
could use a more efficient method which approximates the sum in (14) with
fewer “virtual” degrees as explained in [15], but that is unnecessary here.
Typically we integrate (16) to a stable fixed point and then use pseudo-
arclength continuation to follow the fixed point as parameters are varied,
determining the stability of the fixed point from the eigenvalues of the
linearisation of the dynamics about the fixed point [13, 9]. Periodic
solutions are studied in a similar way by putting a Poincaré section in the
flow (at Re[$b(m+1,t)]=0$) and integrating from this section until the
solution next hits this section. Stability is given by the Floquet multipliers
of the periodic solution.

The usual complex-valued order parameter defined for the network (1) is

$Y(t)=\frac{1}{N}\sum_{j=1}^{N}e^{i\theta_{j}}$ (27)

and for (16) the appropriate measure is

$Z(t)=\sum_{k_{in}}p(k_{in})\bar{b}(k_{in},t)$ (28)

We first vary $\epsilon$ with $\sigma=0.01$. Results are shown in Fig. 2,
where panel (a) shows results from (16). For small $\epsilon$ (16) has a
stable fixed point at which the network is incoherent, with $|Z|$ being small.
As $\epsilon$ is increased the fixed point undergoes a subcritical Hopf
bifurcation, becoming unstable. (In the all-to-all coupled network, this Hopf
bifurcation is supercritical.) The unstable periodic orbit created in this
bifurcation is shown with red crosses and it becomes stable in a saddle-node
bifurcation. Thus there is a small range of $\epsilon$ values for which the
network is bistable. (For periodic orbits, the quantity plotted on the
vertical axis is the average over one period of $|Z(t)|$.)

Panel (b) of Fig. 2 shows $|Y|$ for the network (1) with $\epsilon$
quasistatically increased or decreased, using the final state of the network
at one value of $\epsilon$ as the initial condition for the next value. The
value plotted is the mean over 5000 time units of $|Y(t)|$. The bistability
and hysteresis is clear. Networks of the form used in (1) were created using
the configuration model [newman03]. Both self-connections and multiple
connections between oscillators removed in a random way [16]. Values of
$\omega_{j}$ were then assigned as above, using Lorentzian distributions
centred at target frequencies.

Figure 2. (a): $|Z|$ for fixed point (lines) and periodic solutions (symbols)
of (16). Blue solutions are stable while red are unstable. (b): $|Y|$ for the
network (1). Black corresponds to increasing $\epsilon$ and magenta to
decreasing. Parameters: $m=100,M=400,N=2000,\sigma=0.01,\Delta=0.0005$.

Following the Hopf bifurcation and saddle-node bifurcation of periodic orbits
shown in Fig. 2(a) as both $\sigma$ and $\epsilon$ are varied we obtain Fig.
3. Although the range of values of $\epsilon$ for which the system is bistable
is small, it increases as $\sigma$ is increased. We have shown that even with
very different distributions of in-degrees and intrinsic frequencies,
positively correlating them can induce ES in a directed network of Winfree
oscillators.

Figure 3. Curves of Hopf and saddle-node bifurcations as seen in Fig. 2(a) as
both $\sigma$ and $\epsilon$ are varied. The network is bistable between the
curves. Parameters: $m=100,M=400,\Delta=0.0005$.

The results of our numerical investigation into the effects of varying the
correlation between the $k_{in,j}$ and $\omega_{j}$ ($\rho_{k\omega}$) are
shown in Fig 4 (panels C and D). We varied $\rho_{k\omega}$ between 0.6 and
0.95 using the Monte-Carlo degree-swapping method explained above and for each
network we quasistatically increased or decreased the coupling strength
$\epsilon$ and measured the time-averaged value of $|Y|$. This is shown in Fig
4C where we see results similar to those in Fig. 2 — a small region of
bistability. Panel D of Fig 4 shows the fraction of effective frequencies as
the coupling strength is progressively increased, at $\rho_{k\omega}=0.95$.

The distribution of frequencies and their fractional evolution corresponds to
the mean and standard deviation of the Gaussian distribution, and the highest
concentration of frequencies around the mean emerge dominant. However, the
effects of having a large value of $\rho_{k\omega}$ for Gaussian distributed
intrinsic frequencies is minimal compared with that for power law distributed
frequencies, considered next.

Figure 4. Panels A and C: $|Y|$ as $\epsilon$ is quasistatically varied for
networks with different values of $\rho_{k\omega}$. Panels B and D:
Progression of mean-effective frequency fraction distribution over increasing
coupling strength $\epsilon$, for $\rho_{k\omega}$ = 0.95; note distinct
colour scales for effective frequency evolution due to underlying $\omega$
distributions (histograms, insets Panels B & D), and distinct $\epsilon$
scales reflecting hysteresis over different coupling strengths. Panels A and B
show results when the intrinsic frequencies are power law distributed with an
exponent of $-3.07$ and in- and out-degrees drawn from power law distribution
with same exponent and neutral correlation between the in- and out-degrees.
Panels C and D show results for Gaussian distributed frequencies with mean
$1.0$ and $\sigma=0.01$, while in- and out-degrees are power law distributed
with exponent $-2.66$ and the in-degrees are highly correlated with the out-
degrees ($\rho=0.89$).

##  4\. Power law frequency distribution

We keep the power law distribution of degrees (21) and now consider the case
where the target distribution frequency distribution is also power law
distributed and limited between $c$ and $C$, but with the exponent as a
parameter, i.e.

$p_{\omega_{0}}(\omega_{0})=\begin{cases}a_{\omega}/\omega_{0}^{\gamma+1}&c\leq\omega_{0}\leq
C\\\ 0&\mbox{ otherwise}\end{cases}$ (29)

where $a_{\omega}=\gamma c^{\gamma}C^{\gamma}/(C^{\gamma}-c^{\gamma})$. As
above, having created a network we randomly choose target frequencies from the
distribution (29), then assign the smallest target frequency to the oscillator
with the smallest in-degree, all the way up to the largest target frequency
being associated with oscillator with the largest in-degree. The actual
$\omega_{j}$ are then chosen from a Lorentzian with HWHM equal to $\Delta$
centred at the target frequency, as above. The dependence of $\omega_{0}$ on
$k$ is then

$\omega_{0}(k)=\widehat{p}_{\omega_{0}}^{-1}(\widehat{p}(k))$ (30)

where $\widehat{p}_{\omega_{0}}$ is the cumulative distribution function of
$p_{\omega_{0}}$, i.e

$\widehat{p}_{\omega_{0}}(\omega_{0})=\frac{a_{\omega}}{\gamma}\left(\frac{1}{c^{\gamma}}-\frac{1}{\omega_{0}^{\gamma}}\right)$
(31)

###  4.1. Highly correlated degree and frequency

####  4.1.1. Independent degrees

We first consider the case of independent degrees, as in Sec. 3. Thus we
numerically investigate (16) with (14), but using (30). We choose
$\Delta=0.01$, set $c=1,C=6$, and initially choose $\gamma=2$. As in Sec. 3,
for small $\epsilon$ the system has a stable fixed point and this becomes
unstable through a subcritical Hopf bifurcation as $\epsilon$ is increased.
The results are shown in Fig. 5 where we see the periodic orbit created in the
Hopf bifurcation, giving the same scenario as in Fig. 2 (a). Quasistatically
increasing $\epsilon$ the solution of (16) would jump from a fixed point to a
periodic state with amplitude significantly larger than zero. Decreasing
$\epsilon$, the solution would jump from a finite-amplitude periodic orbit to
a fixed point. (As above, for periodic orbits, the quantity plotted on the
vertical axis is the average over one period of $|Z(t)|$.)

Figure 5. $|Z|$ for fixed point (lines) and periodic (symbols) solutions of
(16). Blue solutions are stable while red are unstable. Parameters:
$c=1,C=6,m=100,M=400,\gamma=2,\Delta=0.01$.

The behaviour of a particular realisation of the discrete network (1) is
slightly different, since a fixed point of (16) corresponds to an incoherent
solution of (1) for which $|Y|$ is not constant, having small fluctuations
about an average value. An example of such dynamics is shown in Fig. 6(a),
with $\epsilon=0.92$. (Other parameters have the same values as in Fig. 5.) We
see that most oscillators are oscillating, but with independent phases.
Similarly, a periodic solution of (16) corresponds to a solution of (1) for
which $|Y|$ is nearly periodic, with the vast majority of the oscillators
having the same average frequency, as shown in Fig. 6(b) ($\epsilon=1$). Some
with high in-degree fire at multiples of this frequency and some are unlocked,
giving a state referred to by [1] as “partially locked hybrid states.” See
also [26] for examples of these dynamics.

Figure 6. $\sin{\theta}$ shown in colour for a simulation of (1). The
oscillators are sorted by their in-degree. (a): $\epsilon=0.92$. (b):
$\epsilon=1$. Other parameters:
$N=5000,c=1,C=6,m=100,M=400,\gamma=2,\Delta=0.01$.

However for $\gamma=1$ the results are very different. As seen in Fig. 7(a)
the branch with the lower value of $|Z|$ still undergoes a Hopf bifurcation.
Numerical integration of (16) near the bifurcation does not show a stable
periodic orbit so this bifurcation seems subcritical, creating an unstable
periodic orbit as $\epsilon$ is decreased. This orbit never becomes stable,
and we hypothesise that it is destroyed in a homoclinic bifurcation with the
“middle” unstable branch. Thus the system has no stable periodic orbits, but
rather a region of bistability between two fixed points with different values
of $|Z|$. Either increasing or decreasing $\epsilon$ the network jumps from
one fixed point to another. Fig. 7(b) shows the mean firing rate $F$ across
the network and we see that the branch with large $|Z|$ has small $F$ and vice
versa. So even if $|Z|$ is large, normally indicating synchronous
oscillations, here it corresponds to a state in which most of the oscillators
are locked at an approximate fixed point, not firing.

Figure 7. (a) $|Z|$ and (b) mean firing rate $F$ for fixed points of (16).
Solid: stable; dashed: unstable. The bifurcation on the lower branch in (a) is
a subcritical Hopf. Parameters: $c=1,C=6,m=100,M=400,\gamma=1,\Delta=0.01$.

Fig. 8 shows (on a logarithmic scale) the expected frequency for oscillators
with in-degree $k_{in}$, given by (19), for three coexisting steady states in
Fig. 7 at $\epsilon=1.7$. The stable solution with highest frequencies has the
lowest value of $|Z|$ and vice versa. Interestingly, for the upper curve the
frequency increases with in-degree $k_{in}$, but for the lower two curves the
maximum frequency does not occur at either extreme of the $k_{in}$ values.

Figure 8. Expected frequency for oscillators with in-degree $k_{in}$, given by
(19), for three coexisting steady states at $\epsilon=1.7$. See Fig. 7. Blue:
stable solutions; red: unstable solution. Other parameters:
$c=1,C=6,m=100,M=400,\gamma=1,\Delta=0.01$.

A third scenario occurs for $\gamma=1.5$, as seen in Fig. 9. The fixed point
that is stable for small $\epsilon$ becomes unstable through a subcritical
Hopf bifurcation as $\epsilon$ is increased as in Fig. 5, but the stable
periodic orbit is destroyed in a SNIC bifurcation which occurs at a slightly
lower value of $\epsilon$ than that at which the Hopf bifurcation occurs. So
if $\epsilon$ is slowly increased the network will jump from one fixed point
to another fixed point. But if $\epsilon$ is then decreased the network will
switch from a fixed point to a stable periodic orbit. This stable orbit is
then destroyed in a saddle-node bifurcation as $\epsilon$ is further decreased
and the network will jump to the original fixed point.

Figure 9. $|Z|$ for fixed point (lines) and periodic (symbols) solutions of
(16). Blue solutions are stable while red are unstable. The stable branch of
periodic orbits terminates at the SNIC bifurcation. Parameters:
$c=1,C=6,m=100,M=400,\gamma=1.5,\Delta=0.01$.

####  4.1.2. Identical degrees

We now consider the case of identical in- and out-degrees. We numerically
investigate (16) with (15), using the power law distribution of degrees (21)
and the power law distribution of frequencies (29). The results for $\gamma=2$
are shown in Fig. 10. The scenario is qualitatively the same as in Fig. 9,
with bistability between either two fixed points or between a fixed point and
a periodic solution.

Figure 10. Identical degrees. $|Z|$ for fixed point (lines) and periodic
(symbols) solutions of (16). Blue solutions are stable while red are unstable.
The stable branch of periodic orbits terminates at the SNIC bifurcation.
Parameters: $c=1,C=6,m=100,M=400,\gamma=2,\Delta=0.01$.

For $\gamma=1$ we obtain Fig. 11, showing yet another scenario. Here there are
no Hopf bifurcations, only two saddle-node bifurcations of fixed points, with
region of bistability between them.

Figure 11. Identical degrees. (a) $|Z|$ and (b) mean firing rate $F$ for fixed
points of (16). Solid: stable; dashed: unstable. Parameters:
$c=1,C=6,m=100,M=400,\gamma=1,\Delta=0.01$.

All of the results reported in this section have been verified using
simulations of the discrete network (results not shown).

####  4.1.3. Varying degree-frequency correlation

We repeated our numerical investigation into the effects of varying degree-
frequency correlations — which we denoted $\rho_{k\omega}$ — generated via the
Monte-Carlo swapping and permutation assembly method with power law
distributed frequencies. The results are shown in Fig. 4 panels A and B.
Increasing $\rho_{k\omega}$ corresponds with the appearance and subsequent
widening of the region of hysteresis, for the power law distributed $\omega$.
The Gaussian distribution (Panels C and D) alternatively demonstrates
hysteresis but with narrowing at the highest level of correlation and is more
sensitive to the network realisation, or final assembly of the connectivity
matrix $A$ given the same degree and $\omega$ sequence. Some realisations
demonstrated hysteresis but others did not with all other elements identical,
and compared to the power law distributed frequencies, we observed less
consistency in emergent hysteresis across the multiple network realisations.
Further, the Gaussian $\omega$ results are clearly and simply noisier, partly
due to the scale of plots showing hysteresis at far lower $\epsilon$ values.
This is also due to the modes of synchronisation for these two $\omega$
distributions: the Gaussian with dynamic synchrony and the power law with
destructive or a quiescent synchrony as observed with locking the system at
fixed points of zero frequency. Note the progression of mean effective
frequency distributions with increasing coupling strength (at the highest
$\rho_{k\omega}=0.95$), where the power law $\omega$ system collapses into
quiescence (Fig. 4 panel B) as opposed to the scenario for Gaussian
distributed $\omega$ (Fig. 4 panel D). Coherence frequencies emerge at levels
dictated by the concentration of $\omega$ values for each distribution: power
law close to zero and Gaussian close to the mean ($\mu=1.0$).

##  5\. Discussion

We have considered directed networks of Winfree oscillators with truncated
inverse power law distributions of both in- and out-degrees. We considered the
case of oscillators having independent in- and out-degrees, and also the case
where they are equal. For independent degrees we examined both Gaussian and
power law distributed intrinsic frequencies, having these frequencies highly
correlated with oscillators’ in-degrees as a result of sorting both sets and
associating quantities of equal rank. We investigated the effects of varying
this degree of correlation. For identical in- and out-degrees we also
considered power law distributed frequencies with a positive correlation
between an oscillator’s in-degree and its intrinsic frequency. We varied both
the width of the Gaussian frequency distribution (when used) and the exponent
in the power law frequency distribution and examined the transitions that
occurred as the strength of coupling within the network was varied.

In all cases shown there was an “explosive” transition as the coupling
strength was varied from one fixed point to another, or between a fixed point
and a periodic solution. An exception occurs with lower degree-frequency
correlations for the power law frequency distributed case. A variety of
scenarios were seen. This is in contrast to the similar and much more widely
studied Kuramoto model, for which transitions are only between an incoherent
fixed point and a partially synchronous periodic solution [6, 8]. The range of
scenarios is a result of the form of a Winfree oscillator: for strong enough
coupling an oscillator will approximately lock to a fixed point, so that even
though the measure of synchrony $|Z|$ may increase as coupling strength is
increased, this does not necessarily correspond to synchronous firing, as can
be seen in Figs. 7 and 11.

Acknowledgements: This work was partially supported by the Marsden Fund
Council from Government funding, managed by Royal Society Te Aparangi, grant
number 17-MAU-054.

## References

  * [1] Joel T Ariaratnam and Steven H Strogatz.  Phase diagram for the winfree model of coupled nonlinear oscillators.  Physical Review Letters, 86(19):4278, 2001.
  * [2] S Boccaletti, JA Almendral, S Guan, I Leyva, Z Liu, I Sendiña-Nadal, Z Wang, and Y Zou.  Explosive transitions in complex networks’ structure and dynamics: Percolation and synchronization.  Physics Reports, 660:1–94, 2016.
  * [3] Christoph Börgers and Nancy Kopell.  Synchronization in networks of excitatory and inhibitory neurons with sparse, random connectivity.  Neural Computation, 15(3):509–538, 2003.
  * [4] Aldo Casaleggio, Michael L Hines, and Michele Migliore.  Computational model of erratic arrhythmias in a cardiac cell network: the role of gap junctions.  PLoS One, 9(6):e100288, 2014.
  * [5] Sarthak Chandra, David Hathcock, Kimberly Crain, Thomas M. Antonsen, Michelle Girvan, and Edward Ott.  Modeling the network dynamics of pulse-coupled neurons.  Chaos, 27(3):033102, 2017.
  * [6] Raissa M. D’Souza, Jesus Gómez-Gardeñes, Jan Nagler, and Alex Arenas.  Explosive phenomena in complex networks.  Advances in Physics, 68(3):123–223, 2019.
  * [7] Rafael Gallego, Ernest Montbrió, and Diego Pazó.  Synchronization scenarios in the winfree model of coupled oscillators.  Physical Review E, 96(4):042208, 2017.
  * [8] Jesús Gómez-Gardeñes, Sergio Gómez, Alex Arenas, and Yamir Moreno.  Explosive synchronization transitions in scale-free networks.  Phys. Rev. Lett., 106:128701, Mar 2011.
  * [9] Willy JF Govaerts.  Numerical methods for bifurcations of dynamical equilibria, volume 66.  Siam, 2000.
  * [10] Seung-Yeal Ha, Jinyeong Park, and Sang Woo Ryoo.  Emergence of phase-locked states for the winfree model in a large coupling regime.  Discrete & Continuous Dynamical Systems, 35(8):3417, 2015.
  * [11] Christian Kuehn and Christian Bick.  A universal route to explosive phenomena.  Science Advances, 7(16):eabe3824, 2021.
  * [12] Y. Kuramoto.  Chemical oscillations, waves, and turbulence.  Springer-Verlag, 1984.
  * [13] C. R. Laing.  Numerical bifurcation theory for high-dimensional neural models.  The Journal of Mathematical Neuroscience, 4(1):1, 2014.
  * [14] Carlo R Laing.  Effects of degree distributions in random networks of type-i neurons.  Physical Review E, 103(5):052305, 2021.
  * [15] Carlo R Laing and Christian Bläsche.  The effects of within-neuron degree correlations in networks of spiking neurons.  Biological Cybernetics, 114:337–347, 2020.
  * [16] Carlo R Laing, Christian Bläsche, and Shawn Means.  Dynamics of structured networks of winfree oscillators.  Front Syst Neurosci, 15:631377, 2021.
  * [17] I Leyva, A Navas, I Sendina-Nadal, JA Almendral, JM Buldú, M Zanin, D Papo, and S Boccaletti.  Explosive transitions to synchronization in networks of phase oscillators.  Scientific reports, 3(1):1–5, 2013.
  * [18] Weiqing Liu, Ye Wu, Jinghua Xiao, and Meng Zhan.  Effects of frequency-degree correlation on synchronization transition in scale-free networks.  EPL (Europhysics Letters), 101(3):38002, 2013.
  * [19] Tanushree B Luke, Ernest Barreto, and Paul So.  Complete classification of the macroscopic behavior of a heterogeneous network of theta neurons.  Neural Computation, 25:3207–3234, 2013.
  * [20] Shawn A Means, Christian Bläsche, and Carlo R Laing.  A permutation method for network assembly.  PLoS One, 15(10):e0240888, 2020.
  * [21] Sung Moon, Katherine Cook, Karthikeyan Rajendran, I. G. Kevrekidis, Jaime Cisternas, and C. R. Laing.  Coarse-grained clustering dynamics of heterogeneously coupled neurons.  Journal of Mathematical Neuroscience, 5(1):2, 2015.
  * [22] Theoden I. Netoff, Matthew I. Banks, Alan D. Dorval, Corey D. Acker, Julie S. Haas, Nancy Kopell, and John A. White.  Synchronization in hybrid neuronal networks of the hippocampal formation.  Journal of Neurophysiology, 93(3):1197–1208, 2005.  PMID: 15525802.
  * [23] M.E.J. Newman.  The structure and function of complex networks.  SIAM Review, 45(2):167–256, 2003.
  * [24] E. Ott and T.M. Antonsen.  Low dimensional behavior of large systems of globally coupled oscillators.  Chaos, 18:037113, 2008.
  * [25] E. Ott and T.M. Antonsen.  Long time evolution of phase oscillator systems.  Chaos, 19:023117, 2009.
  * [26] Diego Pazó and Ernest Montbrió.  Low-dimensional dynamics of populations of pulse-coupled oscillators.  Physical Review X, 4(1):011009, 2014.
  * [27] Diego Pazó, Ernest Montbrió, and Rafael Gallego.  The winfree model with heterogeneous phase-response curves: analytical results.  Journal of Physics A: Mathematical and Theoretical, 52(15):154001, 2019.
  * [28] V. Pernice, B. Staude, S. Cardanobile, and S. Rotter.  How structure determines correlations in neuronal networks.  PLoS Computational Biology, 7(5), 2011.  cited By 117.
  * [29] A. Pikovsky, M. Rosenblum, and J. Kurths.  Synchronization.  Cambridge University Press, 2001.
  * [30] Dietmar Plenz and Stephen T Kitai.  Up and down states in striatal medium spiny neurons simultaneously recorded with spontaneous activity in fast-spiking interneurons studied in cortex–striatum–substantia nigra organotypic cultures.  Journal of Neuroscience, 18(1):266–283, 1998.
  * [31] Juan G Restrepo and Edward Ott.  Mean-field theory of assortative networks of phase oscillators.  EPL (Europhysics Letters), 107(6):60006, 2014.
  * [32] Nathan W Schultheiss, Astrid A Prinz, and Robert J Butera.  Phase response curves in neuroscience: theory, experiment, and analysis.  Springer, 2011.
  * [33] Per Sebastian Skardal and Alex Arenas.  Disorder induces explosive synchronization.  Phys. Rev. E, 89:062811, Jun 2014.
  * [34] Per Sebastian Skardal, Juan G Restrepo, and Edward Ott.  Frequency assortativity can induce chaos in oscillator networks.  Physical Review E, 91(6):060902, 2015.
  * [35] Moira L Steyn-Ross, D Alistair Steyn-Ross, James W Sleigh, and DTJ Liley.  Theoretical electroencephalogram stationary spectrum for a white-noise-driven cortex: evidence for a general anesthetic-induced phase transition.  Physical Review E, 60(6):7299, 1999.
  * [36] S.H. Strogatz.  From Kuramoto to Crawford: exploring the onset of synchronization in populations of coupled oscillators.  Physica D, 143:1–20, 2000.
  * [37] Yasuhiro Tsubo, Jun-nosuke Teramae, and Tomoki Fukai.  Synchronization of excitatory neurons with strongly heterogeneous phase responses.  Physical review letters, 99(22):228101, 2007.
  * [38] CJ Wilson and Y Kawaguchi.  The origins of two-state spontaneous membrane potential fluctuations of neostriatal spiny neurons.  Journal of Neuroscience, 16(7):2397–2410, 1996.
  * [39] Arthur T Winfree.  Biological rhythms and the behavior of populations of coupled oscillators.  Journal of Theoretical Biology, 16(1):15–42, 1967.
  * [40] Jinshan Xu, Shakti N Menon, Rajeev Singh, Nicolas B Garnier, Sitabhra Sinha, and Alain Pumir.  The role of cellular coupling in the spontaneous generation of electrical activity in uterine tissue.  PLoS One, 10(3):e0118443, 2015.
  * [41] Xiyun Zhang, Stefano Boccaletti, Shuguang Guan, and Zonghua Liu.  Explosive synchronization in adaptive and multilayer networks.  Physical review letters, 114(3):038701, 2015.
  * [42] Xiyun Zhang, Yong Zou, Stefano Boccaletti, and Zonghua Liu.  Explosive synchronization as a process of explosive percolation in dynamical phase space.  Scientific reports, 4(1):1–6, 2014.