# LUMINOSITY LEVELLING TECHNIQUES FOR THE LHC

B. Muratori<EMAIL_ADDRESS>STFC Daresbury Laboratory  ASTeC and
Cockcroft Institute  UK
T. Pieloni<EMAIL_ADDRESS>CERN  Geneva  Switzerland

###### Abstract

We present the possibilities for doing luminosity levelling at the LHC. We
explore the merits and drawbacks of each option and briefly discuss the
operational implications. The simplest option is levelling with an offset
between the two beams. Crab cavities may also be used for levelling, as may a
squeezing of the beam. There is also the possibility of using the crossing
angle in order to do luminosity levelling. All of these options are explored,
for the LHC and other possible new projects, together with their benefits and
drawbacks.

##  1 Introduction

One of the main measures of a collider’s performance is its luminosity.
However, from the point of view of experiments, what is most important is not
the peak luminosity but, rather, the integrated luminosity. For the detection
of events, it is also preferable that the luminosity remain constant for as
long as possible. Therefore, luminosity levelling can be introduced. This
means that the natural decay of the luminosity is pre-empted and the
luminosity is spoilt initially with respect to the nominal. Then, as the
luminosity decays, it is spoilt less and less in order that it remain constant
for as long as possible. While doing this, it is still very much worthwhile to
start with as high a luminosity as possible, as this will translate in the
luminosity being constant for a longer amount of time after levelling. To
explain what is meant by this, we consider the expression for the luminosity
in the presence of both an offset and a crossing angle such that the crossing
region is illustrated by Fig. 1 (more details can be found in [1, 2]):

${\cal
L}=\frac{N_{1}N_{2}fN_{b}}{4\pi\sigma_{x}\sigma_{y}}W\mathrm{e}^{\frac{B^{2}}{A}}\frac{1}{\sqrt{1+(\frac{\sigma_{s}}{\sigma_{x}}\tan\frac{\phi}{2})^{2}}},$

with

$A=\frac{\sin^{2}\frac{\phi}{2}}{\sigma_{x}^{2}}+\frac{\cos^{2}\frac{\phi}{2}}{\sigma_{s}^{2}},B=\frac{(d_{2}-d_{1})\sin\frac{\phi}{2}}{2\sigma_{x}^{2}},$

and

$W=\mathrm{e}^{-\frac{1}{4\sigma_{x}^{2}}(d_{2}-d_{1})^{2}}.$

$N_{1}$ and $N_{2}$ are the number of protons per bunch for Beams $1$ and $2$,
respectively; $N_{b}$ is the number of colliding bunches per beam;
$\sigma_{x}$, $\sigma_{y}$, and $\sigma_{s}$ are the transverse and
longitudinal bunch dimensions; $\phi$ is the crossing angle; and $d_{1}$ and
$d_{2}$ are the offsets of Beams $1$ and $2$ with respect to the nominal.

Figure 1: The geometry of the interaction region, with the crossing angle and
the offset of the beams.

Various types of luminosity levelling have been suggested. These are explained
below, together with an analysis of their merits and drawbacks, as well as a
discussion of how they satisfy the requirements, from the point of view of
observation, operations, and the LHC experiments [3, 4, 5, 6]. The main types
of levelling are separation, crab cavities or crossing angle, and $\beta^{*}$
squeeze, as illustrated in Fig. 2.

Figure 2: A sketch of the luminous region for different levelling techniques.

##  2 Levelling with offset

The simplest form of levelling is achieved by introducing an offset between
the two colliding beams. This is straightforward from an operational point of
view, and can be implemented easily and quickly if required. It also makes it
possible to do levelling in all IPs independently, as it is done with a local
orbit bump, and it gives a smaller tune spread, therefore leading to smaller
losses. The average number of p-p collisions per bunch crossing, or pile-up,
that experiments can handle is limited, as different events need to be
distinguished. This is particularly important longitudinally, where the vertex
density is critical and levelling with offset allows for this to be kept
constant.

Luminosity levelling with offset suffers from several drawbacks. The most
obvious one is that a different separation leads to a different beam–beam
force being experienced. Therefore, the effect of one beam–beam encounter with
a separation of a few orders of the r.m.s. beam size can change the tune
spread appreciably, as shown in Fig. 3. This leads to a decrease in the extent
of the stability region as a function of the offset, as shown in Fig. 4.

Figure 3: The tune footprints for different separations. Figure 4: The real
tune shift (the extent of the stability region) as a function of the offset.

Note that the minimum is hardly visible and lies just below the big maximum
shown in white. The position and amplitude of this minimum depends on the
collision schedule, the bunch intensities, emittances, octupole settings, and,
in particular, the transverse offset at all the different IPs and the bunches
experiencing head-on collisions. The fewer the head-on collisions, the smaller
is the stability area shown in Fig. 4. In fact, as can be seen from the
figure, there exists a critical separation at which the instability diagram is
a minimum. Therefore, it is believed that levelling with transverse offset
leads to serious complications in ensuring the stability of all bunches
involved. Operationally [7], it is believed that such effects have already
been observed in IP8 where, during a few fills, bunches that were colliding
only in IP8 were lost and suffered substantial reductions in intensity, as
shown in Fig. 5, together with the full separation inferred from the measured
luminosity in IP1 and IP8, given in Fig. 6. Clearly, it is expected that this
effect will become even worse when bunches collide in more IPs and not all of
them experience head-on collisions. Other drawbacks include the fact that the
tune shift keeps changing as the beams are brought into and out of collision,
and that bunches become more sensitive to instabilities with respect to head-
on collisions. The mere fact of going into collision with a separation could
in itself give rise to instabilities or maximize their effect. Finally, there
is a possible emittance growth resulting from the offsets used, as can be seen
in Fig. 7.

Figure 5: The intensity of the IP8 bunches. Figure 6: Full separation in IP8.
Figure 7: Emittance as a function of beam separation.

##  3 Crab Cavity Levelling

Crab cavities have been used successfully, in electron colliders, to increase
the luminosity back to the nominal value in the presence of a crossing angle.
In a similar way, they may be used to perform luminosity levelling by
‘spoiling’ the luminosity initially, by artificially ‘anti’-crabbing the beam,
and subsequently by correcting for the natural exponential decrease in
luminosity through the usual crabbing of the beam. The main advantages of
using crab cavities are that all IPs are independent, and that it is possible
to go back and forth easily by just changing the voltage of the cavities.

Luminosity levelling with crab cavities suffers from several drawbacks. The
most obvious one is that there is, so far, no experience with applying crab
cavities to proton bunches at all, and this would most likely lead to
additional problems from an operational point of view. The longitudinal vertex
density changes with the levelled angle, giving rise to all the problems that
were discussed above for offset levelling. Further, the tunes change with the
crossing angle, and additional noise could be introduced on the colliding
beams, hence reducing the reachable $\xi_{bb}$. Also, the jitter coming from
the cavities needs to be dealt with. Differential phase jitter causes the two
bunches to have a height mismatch, which can significantly reduce luminosity
or cause the bunches to miss. Phase jitter means that the entry time of the
centre of the bunch to is different for the cavities; hence d$x$ is different
for the two beams, as may be seen in Fig. 8.

Figure 8: The effect of crab cavity jitter at the IP, where $t_{0}$ is the
time at which the bunch enters the cavity.

Phase jitter between the cavities causes the two beams to be displaced in the
$x$-plane, which can reduce the luminosity of the collision or even cause the
bunches to miss each other completely. Further information about cavity phase
jitter for the ILC can be found in [8].

##  4 $\beta^{*}$ Levelling

Another option for doing luminosity levelling is to start with a beam the
cross-section of which is larger than the nominal and then gradually squeeze
it as the luminosity spontaneously reduces exponentially: this is known as
$\beta^{*}$ levelling. Now, the stability of the beam relies on impedance
modes in the machine being Landau damped. This is done via a tune spread in
the presence of beam–beam or other non-linearities in the machine. In the
absence of colliding beams, octupoles are used to ensure the required tune
spread for damping and, thereby, beam stability [9]. Beam–beam effects may be
safely ignored before the $\beta^{*}$ squeeze; however, during the squeeze,
the $\beta$ function grows dramatically in the region near the IP. This
reduces the separation between the two beams even before they are brought into
collision. The stability region before and after the squeeze with two
different octupole settings is shown in Fig. 9.

Figure 9: The horizontal stability diagram for different octupole settings
before and after the $\beta^{*}$ squeeze.

Clearly, it is preferable to use the positive octupole polarity; however, the
strength required means that there are some detrimental effects associated
with this, namely a reduction in dynamic aperture and feed-down effect. This
is avoided if the squeeze is done when the beams are already colliding head
on. This ensures a much larger tune spread and hence Landau damping, giving a
much larger stability diagram, as shown in Fig. 10.

Figure 10: The horizontal stability diagram with and without head-on
collision.

The principle of $\beta^{*}$ levelling is illustrated in Fig. 11, for an
experiment done in 2012 [7] where the beam was slowly squeezed, as a function
of time, and the luminosity increased. Physically, there is no difference
between this and keeping the luminosity constant, as it naturally degrades.
All parameters, such as beam size, tunes, and orbit, should be monitored while
the luminosity is being squeezed. The main advantages are as follows: there is
a constant longitudinal vertex density for the experiments; the tunes do not
change and are constant over the fill; and it is more stable with the largest
area of Landau damping.

Figure 11: The luminosity evolution during fill $2828$.

As the tune spread from head-on beam–beam does not depend on $\beta^{*}$,
leveling with it would allow a constant stability diagram to be maintained
during the procedure, as opposed to what happens when the levelling is done
with just offset. Figure 12 shows a comparison of the measured luminosity
reduction factors when doing $\beta^{*}$ levelling for the experiment
performed at the LHC in 2012 [7], both at CMS and ATLAS, as well as the
expected reduction.

Figure 12: The measured luminosity reduction factor at CMS and ATLAS, compared
to expectations.

The principal drawbacks are due to the orbit. This has to be kept constant
during the squeeze, as the beams must be kept in collision, which means that a
feed-forward on it is required for robustness from an operational point of
view and could require several changes from normal operations.

##  5 Other Levelling Possibilities

Several other possibilities exist, and these are listed briefly below, Some of
them are very new and have not been fully evaluated yet, while others still
require experimental verification and further studies of their viability.

  * •

_Longitudinal cogging_ : This means introducing a time delay of the order of a
couple of RF periods longitudinally, thereby ensuring that there is only a
partial overlap of the two colliding bunches at the IP. So far, only $1$ or
$5$ RF periods have been implemented experimentally. It appears to be a
relatively easy option to implement; however, it means that levelling is done
at all IPs simultaneously and this is very restrictive for the experiments.
Longitudinal cogging also moves the luminous region longitudinally.

  * •

_The large crossing (Piwinski) angle option_ : This is where the levelling is
actually done with a variation of the crossing angle. This option also varies
the length of the luminous region according to the crossing angle.

  * •

_The flat beam option_ : This has been proposed recently [10] and involves
doing the levelling in one plane only, the same as the crossing angle plane.
This means that the tune shift in the other plane can be kept constant and the
collimators do not have to move much, which could otherwise lead to safety
issues.

##  6 Discussion

Various scenarios for luminosity levelling, all valid working options, have
been presented and their merits and drawbacks have been discussed. The easiest
to implement is the offset option; however, this could lead to instabilities.
Crab cavities introduce an additional complexity, which could turn out to be
very non-trivial. There is also no experience of crab cavities and proton
bunches, and the cavities could also introduce a substantial jitter. Several
other options, such as levelling with a crossing angle, flat beam options, or
longitudinal cogging, have also been discussed. However, it appears that
$\beta^{*}$ levelling, possibly together with some offset as well, is the most
promising option, as it appears to satisfy most of the requirements, both from
the experiments’ point of view and operationally. The main problem with
$\beta^{*}$ levelling is that the orbit has to be kept constant as the
levelling is being done, and this may be rather complex from an operational
point of view.

Ultimately, the most important thing that will determine exactly how the
luminosity shall be levelled, both at the LHC and in possible new projects
[11], is the operational simplicity with which the method can be implemented,
together with the experimental requirements and constraints.

##  7 Acknowledgements

The authors would like to thank the Operation group for all the support during
the testing of the $\beta^{*}$ levelling experiments, in particular W. Herr,
J. Wenninger, S. Redaelli, and M. Lamont.

## References

  * [1] W. Herr and B. Muratori, “Concept of Luminosity,” in Proc. CERN Accelerator School, Intermediate Level Course, Zeuthen, Germany, 2003, CERN-2006-002 (2006).
  * [2] B. Muratori, “Luminosity and Luminous Region Calculations for the LHC,” LHC Project Note 301 (2002).
  * [3] G. Papotti, “Observations of Beam–Beam Effects in the LHC,” these proceedings.
  * [4] R. Jacobsson, “Needs and Requirements from the LHC Physics Experiments,” these proceedings.
  * [5] D. Jacquet, “Implementation and Experience with Luminosity Levelling with Offset Beams,” these proceedings.
  * [6] R. Giachino, “Diagnostics Needs for Beam–Beam Studies and Optimization,” these proceedings.
  * [7] X. Buffat et al., “Results of $\beta^{*}$ Luminosity Leveling MD,” CERN-ATS-Note-2012-071 MD (2012).
  * [8] A. Dexter et al., “ILC Crab Cavity Phase Control System Development and Synchronisation Testing in a Vertical Cryostat Facility,” EUROTeV-Report-2008-073 (2008).
  * [9] J.P. Koutchouk and F. Ruggiero, “A Summary on Landau Octupoles for the LHC,” LHC Project Note 163 (1998).
  * [10] A. Burov, “Circular Modes,” these proceedings.
  * [11] T. Pieloni, “Beam–Beam Studies in the LHC and New Projects,” these proceedings.


# RIDS: Robust Identification of Sparse Gene Regulatory Networks from
Perturbation Experiments

Hoi-To Wai, Anna Scaglione, Uzi Harush, Baruch Barzel and Amir Leshem H.-T.
Wai and A. Scaglione are with School of ECEE, Arizona State University, USA.
Emails: {htwai<EMAIL_ADDRESS>U. Harush and B. Barzel are with
Department of Mathematics, Bar-Ilan University, Israel. Emails:
<EMAIL_ADDRESS>A. Leshem is with Faculty of Engineering,
Bar-Ilan University, Israel. Email<EMAIL_ADDRESS>
###### Abstract

Reconstructing the causal network in a complex dynamical system plays a
crucial role in many applications, from sub-cellular biology to economic
systems. Here we focus on inferring gene regulation networks (GRNs) from
perturbation or gene deletion experiments. Despite their scientific merit,
such perturbation experiments are not often used for such inference due to
their costly experimental procedure, requiring significant resources to
complete the measurement of every single experiment. To overcome this
challenge, we develop the Robust IDentification of Sparse networks (RIDS)
method that reconstructs the GRN from a small number of perturbation
experiments. Our method uses the gene expression data observed in each
experiment and translates that into a steady state condition of the system’s
nonlinear interaction dynamics. Applying a sparse optimization criterion, we
are able to extract the parameters of the underlying weighted network, even
from very few experiments. In fact, we demonstrate _analytically_ that, under
certain conditions, the GRN can be perfectly reconstructed using
$K=\Omega(d_{max})$ perturbation experiments, where $d_{max}$ is the maximum
in-degree of the GRN, a small value for realistic sparse networks, indicating
that RIDS can achieve high performance with a scalable number of experiments.
We test our method on both synthetic and experimental data extracted from the
DREAM5 network inference challenge. We show that the RIDS achieves superior
performance compared to the state-of-the-art methods, while requiring as few
as $\sim 60\%$ less experimental data. Moreover, as opposed to almost all
competing methods, RIDS allows us to infer the directionality of the GRN
links, allowing us to infer empirical GRNs, without relying on the commonly
provided list of transcription factors.

##  1 Introduction

Our ability to understand, predict and manipulate the behavior of sub-cellular
systems is limited by our currently partial maps of the cell’s biochemical
process, obstructing the exposure of microscopic biological mechanisms and
hindering drug development [10, 7, 1, 18, 2]. To make advances, researchers
rely on high-throughput profile technologies, such as DNA microarrays, where a
large amount of gene expression data are used for reverse engineering gene
regulatory networks (GRNs). The scale of these experiments allows one to infer
the GRN by extracting the correlations [21] or other related statistical
similarity measures [20, 11] in the expression levels genes under the
different experimental conditions. Alternative methods rely on techniques such
as feature selection, time series regression, etc., [17, 16, 6, 32, 27, 13,
23, 31, 28, 22]. The problem is that such methods do not directly infer the
causal relationships between genes, often falsely identifying indirect
correlations as physical interactions [3]. Another challenge is that
statistical similarity measures are symmetric, and hence lack information on
the directionality of the links. Therefore, one must rely on exogenous
information, such as a predefined list of transcription factors [21, 15],
without which it is not possible to deduce the direction of the inferred
links.

To address these challenges, researchers attempt to infer the structure of the
GRN directly from genetic perturbation data [8, 6, 32, 12, 19, 22, 26]. In
each experiment one or more genes is perturbed, e.g., knocked out or
over/under-expressed, and the change in the expression of all remaining genes
is observed. Using perturbation data has two main advantages: (i) as opposed
to correlations and other related measures, that may result from spurious
statistical dependencies, perturbations represent a controlled experimental
measure of causal effects [19, 5]; (ii) perturbations provide directional
information, as $i$ will only respond to $j$ if there is a directed link/path
leading from $j$ to $i$. The problem is that perturbation experiments are of
limited scale, with the most comprehensive datasets incorporating at most
several hundreds of perturbations, orders of magnitude less than the $\sim
10^{4}$ genes in a typical GRN. Hence, despite their crucial advantages, we
face severe limitations in our ability to use perturbation data for sub-
cellular network inference.

Here we exploit the fact that real GRNs are extremely sparse, with the typical
gene regulated by at most a few other genes. To take advantage of this, we
develop the Robust IDentification of Sparse networks (RIDS) method, that
allows us to infer GRN from a small number of perturbations. Our method takes
as input a set of perturbation experiments, and as output provides a minimal
(sparse) network $\\{{A}_{ij}\\}_{i,j=1}^{n}$ that is consistent with the
observed data. Testing our approach against the DREAM5 gold standard [21] we
find that with as few as $\sim 35\%$ of the total input provided in DREAM5, we
are able to match, and in certain cases, even supass the top performing
inference methods. Since perturbations also encrypt information on the
direction of the causal influence relationship, RIDS allows us to infer the
directionality of the links, as we demonstrate by applying our method without
using the DREAM5 provided list of transcription factors. Strikingly, we find
that even without the transcription factor list, and despite using
significantly less input data, RIDS remains competitive with the top
preforming methods of DREAM5. Finally, an additional advantage that our
perturbation-based inference provides is that the inferred link weights
$A_{ij}$ capture the true strength of the interactions, providing an estimate
for the rate constants for the different genetic reactions. This is in
contrast to the majority of inference methods that only provide likelihood
estimates, that, on their own, bare no inherent meaning, other than indicating
the likelihood for the existence or in-existence of a link.

Notations. For any natural number $n\in\mathbb{N}$, we denote $[n]$ as the set
$\\{1,2,...,n\\}$. Vectors (resp. matrices) are denoted by boldfaced letters
(resp. capital letters). We denote $x_{i}$ as the $i$th element of the vector
${\bm{x}}$, $[{\bm{E}}]_{{\mathcal{}S},:}$ (resp.
$[{\bm{E}}]_{:,{\mathcal{}S}}$) denotes the submatrix of
${\bm{E}}\in\mathbb{R}^{m\times n}$ with only the rows (resp. columns)
selected from ${\mathcal{}S}\subseteq[m]$ (resp. ${\mathcal{}S}\subseteq[n]$).
Vector ${\bf e}_{k}\in\mathbb{R}^{n}$ is a unit vector with zeros everywhere
except for the $k$th coordinate. The superscript $(\cdot)^{\top}$ denotes
matrix/vector transpose. $\|\cdot\|_{2}$ denotes the Euclidean norm and
$\|\cdot\|_{1}$ is the $\ell_{1}$-norm.

Figure 1: Conceptual overview of RIDS. The proposed method takes the _steady
state_ gene expression levels resulted from $K$ ($K\ll n$ where $n$ is the
number of genes) distinct perturbation experiments and tackle a sparse
optimization problem to recover the GRN. For perfect recovery, under some
conditions on the graph structure, the proposed solution requires only
$K=\Omega(d_{max})$ experiments, where $d_{max}$ is the maximum in-degree of
the GRN.

##  2 Method & Main Results

The GRN of interest is modeled as a directed graph $G=(V,E)$, where
$V=[n]=\\{1,...,n\\}$ is the set of genes and $E\subseteq[n]\times[n]$ is the
edge set. Gene $j$ is said to be a regulator (or deregulator) of gene $i$ if
$(i,j)\in E$, i.e., there exists an edge from $j$ to $i$. The GRN is
associated with a weight matrix ${\bm{A}}\in\mathbb{R}^{n\times n}$ that
encodes the strengths of regulation between the genes. We have
$[{\bm{A}}]_{ij}=A_{ij}=0$ if and only if $(i,j)\notin E$. Here we assume the
absence of self-links, hence ${\rm diag}({\bm{A}})={\bm{0}}$. The time
dependent expression level of gene $i$, $x_{i}(t)$, follows the nonlinear
dynamic:

$\textstyle\dot{x}_{i}(t)=-x_{i}(t)+\sum_{j=1}^{n}A_{ij}\cdot
h(x_{j}(t);{\bm{b}})\;,~{}\forall~{}i\in[n]\;.$ (1)

where $\dot{x}_{i}(t)\mathrel{\mathop{:}}=dx_{i}(t)/dt$ is the rate of change
of the expression level $x_{i}(t)$. The first term on the right hand side
captures the gene $i$’s self dynamics, capturing processes such as
degradation, and the sum captures the impact of $i$’s interacting partners;
$h(x;{\bm{b}})$ is the nonlinear continuous response function describing the
regulatory mechanism such that ${\bm{b}}$ describes its parameters. For
instance, setting $h(x)=c\cdot x^{b}$ describes chemical activation, where
according to the law of mass action $b$ is the level of cooperation in the
activating process. Another frequently used response function is the Hill
function $h(x)=x^{b}/(1+x^{b})$, a saturating function such that
$\lim_{x\rightarrow\infty}h(x)=1$, which captures a _switch-like_ activation
process.

In matrix form, we can express any _unperturbed_ steady state
$\overline{\bm{x}}$ of the system (1) as:

$\overline{\bm{x}}={\bm{A}}{\bm{h}}(\overline{\bm{x}};{\bm{b}})\;,\vspace{-.1cm}$
(2)

where
${\bm{h}}({\bm{x}};{\bm{b}})\mathrel{\mathop{:}}=(h(x_{1};{\bm{b}}),~{}h(x_{2};{\bm{b}}),\ldots,~{}h(x_{n};{\bm{b}}))^{\top}$
is a column vector. Next, we consider a set of $K$ distinct perturbation
experiments. In each experiment $k=1,2,...,K$, we _fix_ the state of gene $k$
at a desired value $z_{k}\in\mathbb{R}$, i.e., $x_{k}(t)=z_{k}$ for all $t$.
The perturbed steady state ${\bm{x}}[k]$ takes the form

${\bm{x}}[k]=({\bm{I}}-{\bf e}_{k}{\bf
e}_{k}^{\top}){\bm{A}}{\bm{h}}({\bm{x}}[k];{\bm{b}})+z_{k}{\bf e}_{k}\;.$ (3)

Notice that when $z_{k}=0$, this corresponds to the deletion of gene $k$ in
the experiment.

Suppose that the parameter ${\bm{b}}$ is known. In Eq. (3), the values of all
${\bm{x}}[k]$ are extracted from the experimental results, hence despite the
nonlinear interactions, the equation is linear in the unknown ${\bm{A}}$. As
for each $k$, Eq. (3) constitutes a set of $(n-1)$-_linear_ equations in the
unknown parameter ${\bm{A}}$. Naively, we must conduct perturbation
experiments for all $n$ genes, yielding the required $n(n-1)$ equations to
reconstruct the $n\times(n-1)$ off-diagonal terms of the unknown GRN. Such
comprehensive perturbation experiments, however, are seldom available. Indeed,
the GRN of most organisms comprises $n\sim 10^{3}$ genes, far exceeding the
scale of the majority of microarray experiments, which, given the level of
available resources, consists of $K\sim 10^{1}-10^{2}$ experiments. Hence, we
shall focus on the limit where $K\ll n$ and derive a sparse optimization
algorithm to extract ${\bm{A}}$ from the resulting underdetermined linear
system (3).

We observe a curious property that can help revealing the structure of
${\bm{A}}$. Denote ${\bm{a}}_{col,k}$ as the $k$th column vector of
${\bm{A}}$. Let ${\nabla}{\bm{h}}({\bm{x}}[k]);{\bm{b}})$ be the diagonal
matrix with the $i$th diagonal element being
$h^{\prime}({x}_{i}[k];{\bm{b}})$, i.e., the derivative of
$h(x_{i}[k];{\bm{b}})$ with respect to $x_{i}[k]$ evaluated at $x_{i}[k]$, we
have:

###### Proposition 1

Consider the dynamics (1). Assume that the perturbation in the steady states
for the $k$th experiment, $\overline{\bm{x}}-{\bm{x}}[k]$, is small and
$\lambda_{max}(({\bm{I}}-{\bf e}_{k}{\bf
e}_{k}^{\top}){\bm{A}}{\nabla}{\bm{h}}({\bm{x}}[k];{\bm{b}}))<1$. The
perturbation in the steady states can be approximated by:

$\overline{\bm{x}}-{\bm{x}}[k]\approx([\overline{\bm{x}}]_{k}-z_{k}){\bf
e}_{k}+([\overline{\bm{x}}]_{k}-z_{k})h^{\prime}(z_{k}){\bm{a}}_{col,k}\;.\vspace{-.1cm}$
(4)

The proof can be found in Appendix A. Proposition 1 implies that the
perturbation introduced by the $k$th experiment is _limited_ only to the
_direct out-neighbor_ of the perturbed node $k$. This matches the result in
[4], which showed that the influence of a perturbation on a node decays
exponentially fast with respect to the shortest distance to the perturbed
node. In light of Proposition 1, we let $\delta>0$ be a pre-defined threshold
and consider the index set:

${\mathcal{}S}=\bigcup_{j=1}^{K}\Big{\\{}(i,j)\in[n]\times[n]:\frac{\big{[}\overline{\bm{x}}-{\bm{x}}[j]\big{]}_{i}}{[\overline{\bm{x}}]_{j}-z_{j}}<\delta\Big{\\}}\;.$
(5)

We also define ${\mathcal{}S}_{i}$ as the restriction of ${\mathcal{}S}$ to
the $i$th row of ${\bm{A}}$, notice that ${\mathcal{}S}_{i}\subseteq[K]$.
Importantly, ${\mathcal{}S}$ can be treated as an estimate for the locations
of _zeros/non-edges_ in the GRN ${\bm{A}}$.

Let $(\overline{\bm{x}},\\{{\bm{x}}[k]\\}_{k=1}^{K})$ be the _response matrix_
in which we gather the gene expression data from $K$ perturbation experiments.
When ${\bm{b}}$ is known, we can recover ${\bm{A}}$ by solving for each
$i\in[n]$:

$\displaystyle\displaystyle\min_{\hat{\bm{a}}_{i}}\displaystyle~{}~{}\|\hat{\bm{a}}_{i}\|_{1}~{}~{}{\rm
s.t.}$
$\displaystyle~{}~{}[\overline{\bm{x}}]_{i}=\hat{\bm{a}}_{i}^{\top}{\bm{h}}(\overline{\bm{x}};{\bm{b}}),~{}({\bm{x}}[k])_{i}=\hat{\bm{a}}_{i}^{\top}{\bm{h}}({\bm{x}}[k];{\bm{b}}),~{}\forall~{}k\in[K]\setminus\\{i\\}$
(6a)
$\displaystyle~{}~{}[\hat{\bm{a}}_{i}]_{i}={0},~{}[\hat{\bm{a}}_{i}]_{j}=0,~{}\forall~{}j\in{\mathcal{}S}_{i},~{}\hat{\bm{a}}_{i}\in\mathbb{R}^{n}\;.$
(6b)

The solution to the above problem, $\hat{\bm{a}}_{i}$, serves as an estimate
for the $i$th row of ${\bm{A}}$, i.e., ${\bm{a}}_{i}$. We emphasize that the
sparsity pattern of the optimum $\hat{\bm{a}}_{i}$,
$\\{j\in[n]:[\hat{\bm{a}}_{i}]_{j}\neq 0\\}$, can be interpreted as the set of
regulators/deregulators to gene $i$ and the ranked list of non-zero elements
of $\hat{\bm{a}}_{i}$ indicates the probability of an existing edge in the
GRN. The problem can be solved efficiently in parallel using general software
package such as cvx (Available at http://www.cvxr.com/cvx).

###  2.1 A sufficient condition for _perfect_ GRN Recovery

To understand the fundamental limits of recovering the GRN with (6), we study
the scenario when the steady states $\overline{\bm{x}}$ and ${\bm{x}}[k]$ are
obtained _with no noise_ , i.e., they satisfy the equalities (2) and (3). The
challenge in the analysis is that the _undetermined_ linear system (6a)
depends on the true network ${\bm{A}}$ itself which is a sparse matrix, and
the dynamical system is non-linear. We develop a new sparse recoverability
condition given that ${\bm{b}}$ is known. To describe the result, let us state
the simplifying assumptions below.

###### H 1

(a) Set ${\mathcal{}S}$ is the complement of the support of ${\bm{A}}$; (b)
matrix ${\bm{A}}$ is non-negative; (c) the approximation in Proposition 1 is
exact; (d) ${\bm{h}}({\bm{x}};{\bm{b}})$ admits an exact first order Taylor
approximation; (e) $\overline{x}_{k}-z_{k}\geq 0$ for all $k\in[K]$.

###### Theorem 1

Assume H1. For each $i\in[n]$, if the induced bipartite graph
$G({\mathcal{}S}_{i},[n])$ is an $(\alpha,\delta)$-unbalanced expander graph
with left degree bounded in $[d_{l},d_{u}]$ such that

$(1+(d_{l}/d_{u})\delta)\|{\bm{a}}_{i}\|_{0}\leq\alpha
n,~{}~{}2(d_{l}/d_{u})\delta>\sqrt{5}-1\;,\vspace{-.0cm}$ (7)

then solving (6) with the additional constraint $\hat{\bm{a}}_{i}\geq 0$
yields a unique solution that ${\bm{a}}_{i}=\hat{\bm{a}}_{i}$.

The formal definition of an expander bipartite graph can be found in [14, 29]
or in Definition 1 of the appendix. The proof of Theorem 1 is relegated to
Appendix B. A curious fact about Theorem 1 is that the condition (7) depends
on the graph structure of ${\bm{A}}$ as well as the sparsity of each row
${\bm{a}}_{i}$ of ${\bm{A}}$. This is due to the fact that the linear system
(6a) depend on ${\bm{A}}$ itself.

A possible scenario satisfying the conditions in Theorem 1 requires choosing
the set of perturbed genes $[K]$ such that each gene in $V$ is regulated by a
similar number ($\sim d$) of the genes in $[K]$. As we shall argue in Appendix
B.1, fix $\alpha\in(0,1)$, when $|{\mathcal{}S}_{i}|\geq c(d)\cdot\alpha n$
for some $c(d)>1$ then $G({\mathcal{}S}_{i},V)$ is an
$(\alpha,1-1/(d-1))$-expander with $d_{l}=d-1$, $d_{u}=d$ in high probability.
Applying Theorem 1 shows that perfect recovery is possible if

$K\approx|{\mathcal{}S}_{i}|>c(d)\cdot(1+\delta(d-1)/d)\|{\bm{a}}_{i}\|_{0}\geq
c(d)\cdot((2d-2)/d)\|{\bm{a}}_{i}\|_{0},~{}\forall~{}i\in[n]\;,$ (8)

where $\|{\bm{a}}_{i}\|_{0}$ counts the number of non-zeros in the vector
${\bm{a}}_{i}$. For example, when $d=10$, $\|{\bm{a}}_{i}\|_{0}=0.1n$, we have
$c(d)\cdot((2d-2)/d)\approx 2.91$. Since $\|{\bm{a}}_{i}\|_{0}\leq d_{max}$,
where $d_{max}$ is the maximum in-degree, it shows that $K=\Omega(d_{max})$ is
a sufficient condition for perfect recovery.

Finally, we remark that the results in Theorem 1 gives a condition for the
perfect recovery or identification of ${\bm{A}}$. This is significantly more
challenging than merely _reconstructing_ the GRN as required by DREAM5 [21],
which focuses only on detecting the edges $E$ of the GRN.

Figure 2: Overview of the proposed RIDS method applied on empirical data. In
the pre-processing stage, we first recover the unperturbed steady-state
expression levels $\tilde{\overline{\bm{x}}}$ by analyzing the principal
component of the stacked response matrix $\overline{\bm{X}}_{obs}$. Then, an
orthogonal projection is applied to the perturbed steady-state expression
levels to recover ${\bm{x}}[k]-\overline{\bm{x}}$ for each perturbation
condition. This forms $(K+1)$ vectors where each of them correspond to a
distinct perturbation condition (including no perturbation). Finally, we
tackle the robust sparse network recovery problem (12) for GRN recovery via
the AO procedure (13).

###  2.2 Robust IDentification of Sparse networks (RIDS) method

So far, the theoretical model above assumes that the steady state expression
data are measured noiselessly and the parameter ${\bm{b}}$ in the model
response function $h(x;{\bm{b}})$ is known. In practice, ${\bm{b}}$ is unknown
and maybe different for every gene. We have the _noisy_ expression data:

$\tilde{\overline{\bm{x}}}=\overline{\bm{x}}+\overline{\bm{\epsilon}}~{}~{}~{}\text{and}~{}~{}~{}\tilde{\bm{x}}[k]={\bm{x}}[k]+\bm{\epsilon}[k]\;,\vspace{-.1cm}$
(9)

where the vectors $\overline{\bm{\epsilon}}$, $\bm{\epsilon}[k]$ represent
some additive noise of unknown distribution. Let $i\in[n]$, ${\bm{b}}_{i}$ be
the parameter of gene $i$ and define the matrix/vector:

${\bm{y}}_{i}\mathrel{\mathop{:}}=\left(\begin{array}[]{c}\tilde{x}_{i}[1]\\\
\vdots\\\ \tilde{x}_{i}[K]\\\
\tilde{\overline{x}}_{i}\end{array}\right)~{}~{}\text{and}~{}~{}{\bm{H}}_{i}({\bm{b}}_{i})\mathrel{\mathop{:}}=\left(\begin{array}[]{c}{\bm{h}}(\tilde{{\bm{x}}}[1];{\bm{b}}_{i})^{\top}\\\
\vdots\\\ {\bm{h}}(\tilde{{\bm{x}}}[K];{\bm{b}}_{i})^{\top}\\\
{\bm{h}}(\tilde{\overline{\bm{x}}};{\bm{b}}_{i})^{\top}\end{array}\right)\;.$
(10)

Let $\hat{\bm{b}}_{i}$ be an estimate of ${\bm{b}}_{i}$. Naturally, one would
like to relax the equalities in (6a) and minimize
$\lambda\|\hat{\bm{a}}_{i}\|_{1}+\|{\bm{y}}_{i}-{\bm{H}}_{i}(\hat{\bm{b}}_{i})\hat{\bm{a}}_{i}\|_{2}\;$.
However, we observe that the $k$th element of ${\bm{y}}_{i}$ is expressed as:

$\tilde{x}_{i}[k]={\bm{h}}({\bm{x}}[k];{\bm{b}}_{i})^{\top}{\bm{a}}_{i}+\epsilon_{i}[k]={\bm{h}}(\tilde{\bm{x}}[k];\hat{\bm{b}}_{i})^{\top}{\bm{a}}_{i}+\bm{\delta}[k]^{\top}{\bm{a}}_{i}+\epsilon_{i}[k]\;,\vspace{-.1cm}$
(11)

where
$\bm{\delta}[k]\mathrel{\mathop{:}}={\bm{h}}({\bm{x}}[k];{\bm{b}}_{i})-{\bm{h}}(\tilde{\bm{x}}[k];\hat{\bm{b}}_{i})$
is an unknown vector that scales with the magnitude of $\bm{\epsilon}[k]$. The
difference vector ${\bm{y}}_{i}-{\bm{H}}_{i}({\bm{b}}_{i})\hat{\bm{a}}_{i}$ is
dependent on ${\bm{a}}_{i}$ and can not be modeled as an additive noise.

On the other hand, the unknown parameter ${\bm{b}}_{i}$ lies in a parameter
set ${\mathcal{}B}$. As such, our robust identification of sparse networks
(RIDS) method tackles — for each $i\in[n]$:

$\displaystyle\displaystyle\min_{\hat{\bm{a}}_{i},\hat{\bm{b}}_{i}}$
$\displaystyle\displaystyle~{}~{}J(\hat{\bm{a}}_{i};\hat{\bm{b}}_{i})\mathrel{\mathop{:}}=\|{\bm{y}}_{i}-{\bm{H}}_{i}(\hat{\bm{b}}_{i})\hat{\bm{a}}_{i}\|_{2}+\rho\|\hat{\bm{a}}_{i}\|_{1}+\gamma\|\hat{\bm{a}}_{i}\|_{2}$
(12a) $\displaystyle{\rm s.t.}$
$\displaystyle~{}~{}[\hat{\bm{a}}_{i}]_{i}={0},~{}[\hat{\bm{a}}_{i}]_{j}=0,~{}\forall~{}j\in{\mathcal{}S}_{i},~{}\hat{\bm{a}}_{i}\in\mathbb{R}^{n},~{}\hat{\bm{b}}_{i}\in{\mathcal{}B}\;,\vspace{-.1cm}$
(12b)

where $\rho,\gamma>0$ are fixed regularization parameter. The derivation
details of (12) can be found in Appendix C. This formulation is akin to the
matrix uncertainty (MU) selector in [25] for sparse recovery with uncertainty.
Despite being robust to measurement error, the above formulation
simultaneously solves for the best model parameter that fits with the
expression data.

However, (12) is a _non-convex_ problem due to the multiplicative coupling in
the least square objective function. The problem cannot be solved directly
using off-the-shelf packages. The RIDS method applies an alternating
optimization (AO) approach to get around with the issue, i.e., by running the
following procedure for each $i\in[n]$ — let $L$ be the maximum number of
iterations:

$\begin{split}{\sf for}~{}~{}&\ell=1,2,3,\ldots,L\\\
&\textstyle\hat{\bm{a}}_{i}^{\ell+1}\leftarrow\arg\min_{\hat{\bm{a}}_{i}}~{}J(\hat{\bm{a}}_{i};\hat{\bm{b}}_{i}^{\ell})~{}~{}{\rm
s.t.}~{}~{}\eqref{eq:sense_r_b}~{}{\sf satisfied}\;,\\\
&\textstyle\hat{\bm{b}}_{i}^{\ell+1}\leftarrow\arg\min_{{\bm{b}}}~{}\|(\hat{\bm{b}}_{i}^{\ell}-\epsilon\cdot{\nabla}_{\bm{b}}J(\hat{\bm{a}}_{i}^{\ell+1};{\bm{b}}^{\ell}))-{\bm{b}}\|_{2}~{}~{}{\rm
s.t.}~{}~{}{\bm{b}}\in{\mathcal{}B}\;,\end{split}$ (13)

where $\epsilon>0$ is a fixed step size and
${\nabla}_{\bm{b}}J(\hat{\bm{a}}_{i}^{\ell+1};\hat{\bm{b}}_{i}^{\ell})$ is the
gradient of the cost function. The last step is a _projected gradient_ update
step for $\hat{\bm{b}}_{i}$. The RIDS method is summarized in Figure 2. In the
first stage, we apply a pre-process method to de-noise the experimental data
(see Appendix D); in the second stage, we tackle the robust GRN recovery
problem (12) using the AO procedure in (13).

##  3 Numerical Results & Discussions

Figure 3: Numerical Experiments on synthetic data of networks with $n=100$
nodes, averaged over $100$ Monte-Carlo trials. Two random graph models are
considered — Erdos-Renyi graph with connectivity $p=0.1$ and Random regular
graph with constant degree $d=10$. The model response function used is
$h(x)=x^{0.5}/(1+x^{0.5})$ and the parameters are assumed to be known. The top
figures consider the noiseless observation case. (Top-Left) Area under ROC
(AUROC) against the number of perturbation experiments $K$. (Top-Right) Mean
square error of the recovered $\hat{\bm{A}}$ using the proposed method. Our
proposed RIDS method (via solving (6)) outperforms GENIE3 and TIGRESS over all
range of $K$ and achieves perfect recovery at $K\geq 32$ for the random
regular graph model (as indicated by the MSE plot). Notice that the TIGRESS
method has encountered numerical issues for the case with ER graphs. Perfect
recovery was also observed for the ER model for $\sim 70\%$ of the instances
at $K\geq 32$, but however, as the ER graphs tend to have hubs with high
degree, its average MSE performance will be affected. The bottom figures
consider the noisy observation case. (Bottom-Left) Random Erdos-Renyi graphs
with connectivity $0.1$. (Bottom-Right) Random Regular graphs with constant
degree $d=10$. Again, we observe that the proposed method has outperformed
GENIE3 and TIGRESS for the considered range of $K$.

###  3.1 Synthetic Data

_In silico_ data. We test the models when the GRN $G=(V,E)$ with $n=100$. The
weight matrix ${\bm{A}}$ has entries that are uniformly distributed in
$(0,1]$. We evaluate the steady-state gene expression levels subject to gene
deletion using the $4$th order Runge-Katta method. We include the GENIE3
method [17] and TIGRESS method [16] implemented by their respective authors
for benchmarks. The latter two were the best performing methods in the DREAM5
challenge. All algorithms tested are implemented on MATLAB 2016a. For the
zeros index set ${\mathcal{}S}$ in (5), we set $\delta=0.02$.

Analysis. Figure 3 (Top) compares the mean square error of the recovered
$\hat{\bm{A}}$ and the area under an ROC curve (AUROC) for the recovered
network $G$ versus the number of perturbation experiments $K$. We assume
noiseless measurements in this case and solve (6) to recover the network. We
observe that the proposed method achieves an AUROC of $\geq 0.9$ with $K\geq
14$ perturbation experiments, significantly outperforming the GENIE3 and
TIGRESS methods under similar conditions. Moreover, we see that the proposed
RIDS method has a better performance when the underlying graph is a regular
graph. In particular, with $K=32$ perturbation experiments we have perfect
recovery of both the inferred links and interaction strengths, $A_{ij}$, for
the regular graph model. The above observation coincides with Theorem 1 which
predicts with $2.91\times\max_{i\in[n]}\|{\bm{a}}_{i}\|_{0}$ perturbation
experiments, one can perfectly recover the GRN with our proposed approach when
$n\rightarrow\infty$. Nevertheless, having $K\approx 14$ experiments ($\sim
15\%$ of the total number of genes) is sufficient to yield a good GRN
reconstruction performance.

Figure 3 (Bottom) considers the _noisy_ measurement scenario. With reference
to (9), the elements of $\overline{\bm{\epsilon}},\bm{\epsilon}[k]$ are
independently extracted from normal distributions ${\mathcal{}N}(0,0.1)$ and
${\mathcal{}N}(0,0.01)$. We apply the robust formulation (12) with the
regularizing parameters set to $\rho=10^{-5},\gamma=0.5$, $\gamma=0.05$ for
the case with ${\mathcal{}N}(0,0.1)$ noise and ${\mathcal{}N}(0,0.01)$ noise,
respectively, to recover the network, notice that in this scenario ${\bm{b}}$
is known and problem (12) is solved directly. Comparing the average AUROC
performance shows that our method has a consistently better performance than
GENIE3 and TIGRESS. However, as the noise power grows, the advantage of
applying our method declines.

###  3.2 Empirical Data

_In vivo_ data. To test our methodology against empirical data, we focus on
the reconstruction performance of Escherichia coli (E. coli) and Saccharomyces
cerevisiae (S. cerevisiae) from gene perturbation experiments, using the
highly curated datasets collected for the DREAM5 Network Inference
Challenge111Available at
``https://www.synapse.org/#Synapse:syn2787209/wiki/70351—.. The _E. coli_
(resp.  _S. cerevisiae_) dataset describes $805$ (resp. $536$) vectors of
expression level of $n=4511$ (resp. $n=5950$) anonymized genes under different
experimental conditions; the dataset also lists a subset of $TF=334$ (resp.
$TF=333$) genes that are recognized as known _transcription factors_ (TFs),
some of whom are decoys, namely wrongly labeled as such.

Our method relies on the expression levels in the steady state only.
Therefore, we use only $326$ out of $805$ (resp. $238$ out of $536$) vectors
of expression levels in the _E. coli_ (resp.  _S. cerevisiae_) dataset, i.e.,
about $\sim 40\%$ of the data. The remaining vectors correspond to transient
states, which our method is not designed to treat. Upon grouping and denoising
the dataset using the pre-processing method in Appendix D, we are left with
$K=56$ (resp. $K=7$) vectors of expression levels for _E. coli_ (resp.  _S.
cerevisiae_), each with a distinct perturbation conditions. From these
vectors, we must reconstruct $\bm{A}$, capturing the GRN between the
transcription factors and the genes, a total of $334\times 4511$ (resp.
$333\times 5950$) potential links for _E. coli_ (resp.  _S. cerevisiae_).

We use the following model response function for empirical data:

$h(x;{\bm{b}})=0.75\cdot
x^{b_{2}}/(1+b_{1}x^{b_{2}}),~{}{\mathcal{}B}=\\{{\bm{b}}\in\mathbb{R}^{2}~{}:~{}b_{1},b_{2}\geq
0\\}\;.$ (14)

The above encompasses several kinetic interaction models. When
$b_{1}\rightarrow 0$, the interaction model tends to be that of a chemical
activation process; otherwise the interaction model has a saturating effect
close to a switch-like process; $b_{2}$ controls the saturation rate of the
gene interaction. We apply the RIDS method to learn simultaneously the GRN and
the above model parameters. Further details of the experimental procedure can
be found in Appendix D.

| _E. coli_ | _S. cerevisiae_
---|---|---
Methods |  AUROC |  AUPR |  Score |  AUROC |  AUPR |  Score
TIGRESS [16] |  0.595 |  0.069 |  4.41 |  0.517 |  0.02 |  1.082
GENIE3 [17] |  0.617 |  0.093 |  14.79 |  0.518 |  0.021 |  1.387
RankSum |  0.65 |  0.09 |  24.90 |  0.528 |  0.022 |  6.236
bLARS [27] |  N/A |  N/A |  5.841 |  N/A |  N/A |  7.479
_RIDS_ |  0.6808 |  0.0504 |  32.39 |  0.525 |  0.022 |  4.694
_(top 100k)_ |  $1.69\times 10^{-64}$ |  $9.9\times 10^{-2}$ |  |  $3.84\times 10^{-8}$ |  $2.3\times 10^{-2}$ |
_RIDS (opt. ${\bm{b}},$_ |  0.6823 |  0.0508 |  33.29 |  0.525 |  0.021 |  4.161
_top 100k)_ |  $3.13\times 10^{-66}$ |  $8.5\times 10^{-2}$ |  |  $6.47\times 10^{-8}$ |  $2.9\times 10^{-2}$ |
_RIDS (no TF,_ |  0.6745 |  0.0540 |  29.12 |  0.524 |  0.0221 |  4.298
_top 100k)_ |  $2.78\times 10^{-57}$ |  $2.0\times 10^{-2}$ |  |  $3.70\times 10^{-7}$ |  $4.8\times 10^{-2}$ |
iRafNet [23] |  0.641 |  0.112 |  29.26 |  N/A |  N/A |  N/A
GENIMS [31] |  0.705 |  0.052 |  48.33 |  0.533 |  0.02 |  8.454
_RIDS_ |  0.7573 |  0.0574 |  93.28 |  0.5734 |  0.0252 |  62.64
_(top 500k)_ |  $1.04\times 10^{-184}$ |  $2.7\times 10^{-3}$ |  |  $1.5\times 10^{-119}$ |  $2.38\times 10^{-7}$ |

⋆All values/scores are calculated with the top 100k predictions. Exceptions
are the iRafNet, GENIMS and RIDS (top 500k) in the last three rows, that are
based on the top 200k, all, top 500k predictions, respectively.

Table 1: GRN recovery result on the two _in vivo_ dataset. Scores for RankSum,
GENIE3 and TIGRESS are taken directly from the DREAM5’s leaderboard. Notice
that RankSum is the community integrated prediction. For the RIDS method, the
lower scores are the $p$-value for the AUROC/AUPR metrics. The RIDS method in
the $5$th row uses the median optimized parameters learnt for the two
networks, i.e., $b_{1}=0.047,b_{2}=0.5893$ for _E. coli_ and
$b_{1}=0.5571,b_{2}=0.3749$ for _S. cerevisiae_ and we solve (12) with the
fixed parameters for all genes; the $6$th row corresponds to tackling (12)
_without_ using the list of transcription factors. Our method obtains the best
prediction score for _E. coli_ , and the score for _S. cerevisiae_ is
comparable to state-of-the-art.

Analysis. Table 1 shows the GRN recovery result for the truncated top
_100,000_ and _500,000_ predicted links in the GRN, compared to the Gold
Standard. For each method, we evaluate the AUROC, AUPR and the prediction
score, defined similarly as in [21] by ${\sf
Score}\mathrel{\mathop{:}}=(\log_{10}(p_{AUROC})+\log_{10}(p_{AUPR}))/2$,
where $p_{AUROC}$ and $p_{AUPR}$ are the $p$-values for AUROC/AUPR. The
numerical experiment demonstrates that our robust GRN recovery approach is
able to infer GRN from _few steady-state data_ , while achieving superior
performance to the state-of-the-art. In particular, for the _E. coli_ network,
the RIDS method is the top performer among the compared methods, beating even
the community integrated prediction (RankSum) while using $\sim 60\%$ less
experimental data. Our method gives good performance even the list of TFs are
not provided when tackling (12).

The model parameters learnt using RIDS are plotted in Figure 4. We observe
that the parameters learnt are clustered around $b_{1}\sim 0.05$, $b_{2}\sim
0.6$ for the _E. coli_ network and $b_{1}\sim 0.4$, $b_{2}\sim 0.55$ for the
_S. cerevisiae_ network. Suggesting that the interaction model for the former
is close to a chemical activation process, while the latter is close to being
a combination. The parameters learnt are also similar across different genes,
which seems to indiciate that the parameter space may be reduced by using a
single set of parameters for all genes. This is corroborated by the AUROC/AUPR
performance obtained when the RIDS method is re-applied with the same model
parameter for all genes in Table 1. Additional results on the empirical data
can be found in Appendix E, where we further demonstrate that RIDS can infer
directionality from empirical GRNs.

Figure 4: Model parameters learnt for different networks. (Left) _E. coli_
network. (Right) _S. cerevisiae_ network. The top figures show the scatter
plot of the parameters. The bottom figures show the histograms of parameters
learnt, where the  blue patch is the saturation rate parameter $b_{2}$ and the
red patch is parameter $b_{1}$ which controls if the interaction model is more
of the chemical activation type or the switching type, cf. (14). Meanwhile,
the _S. cerevisiae_ shows a higher variance for the $b_{1},b_{2}$ values
learnt for each gene, this is due to the fact that the in vivo data available
for this network is more scarce than the one for _E. coli_.

Computational Complexity. The problem (6) and the AO procedure in (13) can be
solved by applying generic optimization packages such as cvx. It takes about
25 seconds to recover the synthetic networks in Section 3.1 with $n=100$ genes
on a quad-core labtop computer and about $1100$ seconds for the $n=4511$ genes
_E. Coli._ network in Section 3.2 on a 60-cores server. Further speed up will
be possible if we apply specialized optimization tools.

Conclusions. This work proposes the RIDS method for GRN inference using
perturbation experiments. The RIDS method is developed through modelling the
gene expression data as the outcome of perturbing a nonlinear dynamic system.
To improve robustness, the method first applies a subspace projection method
for denoising the gene expression data, then a sparse and robust estimator is
applied to recover the network. The model parameters of the dynamic system
will also be inferred simultaneously. Our theoretical analysis, conducted
under the assumption that the dynamic’s parameters are known, shows that it is
possible to recover the GRN even when there are only a few sets of
perturbation experiment data available. This is in contrary to the common
belief that it requires a large number of experiments to apply similar
methods. Moreover, our experiments on empirical data shows that the RIDS
method compares favorably to the state-of-the-art methods while requiring
$\sim 60\%$ less data. The RIDS method paves the way to study the dynamics of
gene interactions by having the ability to infer the model parameters, which
is important as studied by [24]. For example, our preliminary result suggests
that the genes in the _E. coli_ and _S. cerevisiae_ networks tend to have a
different interaction model with its neighboring genes.

## References

  * [1] R. Albert.  Scale-free networks in cell biology.  Journal of Cell Science, 118(21):4947, 10 2005.
  * [2] A.-L. Barabasi and Z. N. Oltvai.  Network biology: understanding the cell’s functional organization.  Nat Rev Genet, 5(2):101–113, 02 2004.
  * [3] B. Barzel and A.-L. Barabasi.  Network link prediction by global silencing of indirect correlations.  Nat Biotech, 31(8):720–725, 08 2013.
  * [4] B. Barzel and A.-L. Barabasi.  Universality in network dynamics.  Nat Physics, 9:673–681, October 2013.
  * [5] B. Barzel and O. Biham.  Quantifying the connectivity of a network: The network correlation function method.  Physical Review E, 80(4):046104–, 10 2009.
  * [6] R. Bonneau, D. J. Reiss, P. Shannon, M. Facciotti, L. Hood, N. S. Baliga, and V. Thorsson.  The inferelator: an algorithm for learning parsimonious regulatory networks from systems-biology data sets de novo.  Genome Biology, 7(5):R36, 2006.
  * [7] M. Buchanan.  Networks in Cell Biology.  Cambridge University Press, 2010.
  * [8] X. Cai, J. A. Bazerque, and G. B. Giannakis.  Inference of gene regulatory networks with sparse structural equation models exploiting genetic perturbations.  PLoS Comput Biol, 9(5):e1003068–, 05 2013.
  * [9] E. Candes and T. Tao.  Decoding by Linear Programming.  IEEE Trans. Inf. Theory, 51(12):4203–4215, 2005.
  * [10] R. De Smet and K. Marchal.  Advantages and limitations of current network inference methods.  Nat Rev Micro, 8(10):717–729, 10 2010.
  * [11] N. Friedman, M. Linial, I. Nachman, and D. Pe’er.  Using bayesian networks to analyze expression data.  Journal of computational biology, 7(3-4):601–620, 2000.
  * [12] T. S. Gardner, D. di Bernardo, D. Lorenz, and J. J. Collins.  Inferring genetic networks and identifying compound mode of action via expression profiling.  Science, 301(5629):102, 07 2003.
  * [13] M. Ghanbari, J. Lasserre, and M. Vingron.  Reconstruction of gene networks using prior knowledge.  BMC Systems Biology, 9(1):84, 2015.
  * [14] A. Gilbert and P. Indyk.  Sparse recovery using sparse matrices.  Proceedings of the IEEE, 98(6):937–947, 2010.
  * [15] C. T. Harbison, D. B. Gordon, T. I. Lee, N. J. Rinaldi, K. D. Macisaac, T. W. Danford, N. M. Hannett, J.-B. Tagne, D. B. Reynolds, J. Yoo, E. G. Jennings, J. Zeitlinger, D. K. Pokholok, M. Kellis, P. A. Rolfe, K. T. Takusagawa, E. S. Lander, D. K. Gifford, E. Fraenkel, and R. A. Young.  Transcriptional regulatory code of a eukaryotic genome.  Nature, 431(7004):99–104, 09 2004.
  * [16] A.-C. Haury, F. Mordelet, P. Vera-Licona, and J.-P. Vert.  Tigress: Trustful inference of gene regulation using stability selection.  BMC Systems Biology, 6(1):1–17, 2012.
  * [17] V. A. Huynh-Thu, A. Irrthum, L. Wehenkel, and P. Geurts.  Inferring regulatory networks from expression data using tree-based methods.  PLoS ONE, 5(9), Sept 2010.
  * [18] T. Ideker and R. Sharan.  Protein networks in disease.  Genome Research, 18(4):644–652, 04 2008.
  * [19] T. Kang, R. Moore, Y. Li, E. Sontag, and L. Bleris.  Discriminating direct and indirect connectivities in biological networks.  Proceedings of the National Academy of Sciences, 112(41):12893–12898, 10 2015.
  * [20] R. Küffner, T. Petri, P. Tavakkolkhah, L. Windhager, and R. Zimmer.  Inferring gene regulatory networks by anova.  Bioinformatics, 28(10):1376–1382, 2012.
  * [21] D. Marbach, J. C. Costello, R. Kuffner, N. M. Vega, R. J. Prill, D. M. Camacho, K. R. Allison, M. Kellis, J. J. Collins, and G. Stolovitzky.  Wisdom of crowds for robust gene network inference.  Nat Meth, 9(8):796–804, 08 2012.
  * [22] N. Omranian, J. M. O. Eloundou-Mbebi, B. Mueller-Roeber, and Z. Nikoloski.  Gene regulatory network inference using fused lasso on multiple data sets.  Scientific Reports, 6:20533 EP –, 02 2016.
  * [23] F. Petralia, P. Wang, J. Yang, and Z. Tu.  Integrative random forest for gene regulatory network inference.  Bioinformatics, 31(12):i197–i205, 06 2015.
  * [24] M. Ronen, R. Rosenberg, B. I. Shraiman, and U. Alon.  Assigning numbers to the arrows: Parameterizing a gene regulation network by using accurate expression kinetics.  Proceedings of the National Academy of Sciences, 99(16):10555–10560, 08 2002.
  * [25] M. Rosenbaum and A. Tsybakov.  Sparse recovery under matrix uncertainty.  Annals of Statistics, 38(5):2620–2651, 2010.
  * [26] A. Shojaie, A. Jauhiainen, M. Kallitsis, and G. Michailidis.  Inferring regulatory networks by combining perturbation screens and steady state gene expression profiles.  PLoS ONE, 9(2):e82393–, 02 2014.
  * [27] N. Singh and M. Vidyasagar.  blars: An algorithm to infer gene regulatory networks.  IEEE/ACM Transactions on Computational Biology and Bioinformatics, 13(2):301–314, March 2016.
  * [28] H. M. Tran and S. T. Bukkapatnam.  Inferring sparse networks for noisy transient processes.  Scientific Reports, 6:21963, 2016.
  * [29] H.-T. Wai, A. Scaglione, and A. Leshem.  Active Sensing of Social Networks.  IEEE Trans. Sig. and Inf. Proc. over Networks, 2(3):406–419, 2016\.
  * [30] M. Wang, W. Xu, and A. Tang.  A unique ”nonnegative” solution to an underdetermined system: From vectors to matrices.  IEEE Trans. Signal Process., 59(3):1007–1016, 2011.
  * [31] J. Wu, X. Zhao, Z. Lin, and Z. Shao.  Large scale gene regulatory network inference with a multi-level strategy.  Molecular BioSystems, 12(2):588–597, 2016.
  * [32] K. Y. Yip, R. P. Alexander, K.-K. Yan, and M. Gerstein.  Improved reconstruction of in silico gene regulatory networks by integrating knockout and perturbation data.  PLoS ONE, 5(1):e8121–, 01 2010.

##  Appendix A Proof of Proposition 1

Define

$\overline{\bm{x}}-{\bm{x}}[k]=([\overline{\bm{x}}]_{k}-z_{k}){\bf
e}_{k}+\bm{\epsilon},$ (15)

such that the $k$th component of $\bm{\epsilon}$ is zero. Our goal is to find
$\bm{\epsilon}$. We have:

$\overline{\bm{x}}-{\bm{x}}[k]={\bm{A}}({\bm{h}}(\overline{\bm{x}})-{\bm{h}}({\bm{x}}[k]))+{\bf
e}_{k}{\bf e}_{k}^{\top}{\bm{A}}{\bm{h}}({\bm{x}}[k])+z_{k}{\bf e}_{k}.$ (16)

Notice that ${\bf e}_{k}$ is in the null space of $({\bm{I}}-{\bf e}_{k}{\bf
e}_{k}^{\top})$. Left-multiplying $({\bm{I}}-{\bf e}_{k}{\bf e}_{k}^{\top})$
to the both sides of the equation yields:

$\begin{split}\bm{\epsilon}&=({\bm{I}}-{\bf e}_{k}{\bf
e}_{k}^{\top}){\bm{A}}({\bm{h}}(\overline{\bm{x}})-{\bm{h}}({\bm{x}}[k]))\\\
&\approx({\bm{I}}-{\bf e}_{k}{\bf
e}_{k}^{\top}){\bm{A}}{\nabla}{\bm{h}}({\bm{x}}[k])(([\overline{\bm{x}}]_{k}-z_{k}){\bf
e}_{k}+\bm{\epsilon})\end{split}$ (17)

where we have taken the Taylor’s expansion for ${\bm{h}}(\overline{\bm{x}})$
centered at ${\bm{x}}[k]$ and the approximation is accurate when the
perturbation $\overline{\bm{x}}-{\bm{x}}[k]$ is small. This gives

$\bm{\epsilon}\approx({\bm{I}}-{\bf e}_{k}{\bf
e}_{k}^{\top}){\bm{A}}{\nabla}{\bm{h}}({\bm{x}}[k])\bm{\epsilon}+([\overline{\bm{x}}]_{k}-z_{k})({\bm{I}}-{\bf
e}_{k}{\bf e}_{k}^{\top}){\bm{A}}{\nabla}{\bm{h}}({\bm{x}}[k]){\bf e}_{k}$
(18)

Moreover, we notice that ${\nabla}{\bm{h}}({\bm{x}}[k])$ is a diagonal matrix
with $h^{\prime}(z_{k})$ on its $k$th entry, therefore the latter term can be
simplified as

$\begin{split}([\overline{\bm{x}}]_{k}-z_{k})({\bm{I}}-{\bf e}_{k}{\bf
e}_{k}^{\top}){\bm{A}}{\nabla}{\bm{h}}({\bm{x}}^{\prime}){\bf
e}_{k}&=([\overline{\bm{x}}]_{k}-z_{k})h^{\prime}(z_{k})({\bm{I}}-{\bf
e}_{k}{\bf e}_{k}^{\top}){\bm{a}}_{col,k}\\\
&=([\overline{\bm{x}}]_{k}-z_{k})h^{\prime}(z_{k}){\bm{a}}_{col,k}\;,\end{split}$
(19)

where the last equality is due to the fact that the $k$th element of
${\bm{a}}_{col,k}$ is zero.

Finally, we observe that

$\begin{split}\bm{\epsilon}&=\big{(}{\bm{I}}-({\bm{I}}-{\bf e}_{k}{\bf
e}_{k}^{\top}){\bm{A}}{\nabla}{\bm{h}}({\bm{x}}[k])\big{)}^{-1}\cdot([\overline{\bm{x}}]_{k}-z_{k})h^{\prime}(z_{k}){\bm{a}}_{col,k}\\\
&=\big{(}{\bm{I}}+({\bm{I}}-{\bf e}_{k}{\bf
e}_{k}^{\top}){\bm{A}}{\nabla}{\bm{h}}({\bm{x}}^{\prime})\\\ &\hskip
56.9055pt+(({\bm{I}}-{\bf e}_{k}{\bf
e}_{k}^{\top}){\bm{A}}{\nabla}{\bm{h}}({\bm{x}}^{\prime}))^{2}+\cdots\big{)}\cdot([\overline{\bm{x}}]_{k}-z_{k})h^{\prime}(z_{k}){\bm{a}}_{col,k}\\\
&\approx([\overline{\bm{x}}]_{k}-z_{k})h^{\prime}(z_{k}){\bm{a}}_{col,k}\;,\end{split}$
(20)

where the second equality is due to the Taylor’s expansion. Notice that the
series expansion holds when $\lambda_{max}(({\bm{I}}-{\bf e}_{k}{\bf
e}_{k}^{\top}){\bm{A}}{\nabla}{\bm{h}}({\bm{x}}^{\prime}))<1$.

##  Appendix B Proof of Theorem 1

Proof Outline. Using the assumptions stated in the Theorem, we first reduce
the linear system into a simple form that involves an underdetermined system
with a _sparse_ sensing matrix. The sensing matrix is then found to have the
same structure/support as the bipartite graph formed by the edges _from_ the
perturbed nodes _to_ all other nodes. Finally, the desired perfect recovery
condition is given as a consequence of the expander graph property of this
bipartite graph.

Let us give a formal definition of expander graph:

###### Definition 1

An $(\alpha,\delta)$-unbalanced expander graph, $G(A,B)$, is an
$A,B$-bipartite graph (bigraph) with $|A|=n,|B|=m$ with left degree bounded in
$[d_{l},d_{u}]$, i.e., $d(v_{i})\in[d_{l},d_{u}]$ for all $v_{i}\in A$, such
that for any $S\subseteq A$ with $|S|\leq\alpha n$, we have
$\delta|E(S,B)|\leq|N(S)|$, where $E(S,B)$ is the set of edges from $S$ to $B$
and $N(S)=\\{v_{j}\in B:\exists~{}v_{i}\in S~{}s.t.~{}(v_{j},v_{i})\in E\\}$
is the neighbor set of $S$ in $B$.

We are ready to begin our proof. Let ${\bm{a}}_{i}$ be the $i$th row of
${\bm{A}}$ and define a set of vectors $\\{{\bm{y}}_{i}\\}_{i=1}^{n}$ such
that $[{\bm{y}}_{i}]_{k}\mathrel{\mathop{:}}=\overline{x}_{i}-x_{i}[k]$. We
observe that ${\bm{y}}_{i}$ is a collection of the data points that depend on
${\bm{a}}_{i}$. For simplicity, let us focus on the case when $i\notin[K]$,

$[{\bm{y}}_{i}]_{k}={\bm{a}}_{i}^{\top}({\bm{h}}(\overline{\bm{x}})-{\bm{h}}({\bm{x}}[k]))\approx(\overline{x}_{k}-z_{k})\cdot{\bm{a}}_{i}^{\top}{\nabla}{\bm{h}}(\overline{\bm{x}})\big{(}h^{\prime}(z_{k}){\bm{a}}_{col,k}+{\bf
e}_{k}\big{)}\;,$ (21)

where we have applied Proposition 1 and the first order Taylor approximation
on ${\bm{h}}(\overline{\bm{x}})-{\bm{h}}({\bm{x}}[k])$ to yield the result
above. Notice that we have dropped the dependence on ${\bm{b}}$ as the ODE
parameter is assumed to be known in this setting. After some manipulations and
applying the assumptions stated in the theorem, we can express the equation
above as

${\bm{y}}_{i}=\left(\begin{array}[]{cc}\bm{\Lambda}&{\bm{0}}_{K\times(n-K)}\end{array}\right){\bm{a}}_{i}+\underbrace{\left(\begin{array}[]{c}(\overline{x}_{1}-z_{1})h^{\prime}(z_{1})\cdot{\nabla}{\bm{h}}(\overline{\bm{x}}){\bm{a}}_{col,1}^{\top}\\\
\vdots\\\
(\overline{x}_{K}-z_{K})h^{\prime}(z_{K})\cdot{\nabla}{\bm{h}}(\overline{\bm{x}}){\bm{a}}_{col,K}^{\top}\\\
\end{array}\right)}_{\mathrel{\mathop{:}}=\tilde{\bm{E}}_{i}}{\bm{a}}_{i}\;,$
(22)

where $\bm{\Lambda}$ is an $K\times K$ diagonal matrix with the $k$th element
being $[\bm{\Lambda}]_{kk}=\overline{x}_{k}-z_{k}$. The challenge is that as
the non-zero elements of $\bm{\Lambda}$ has a larger magnitude than the matrix
$\tilde{\bm{E}}_{i}$ in the latter matrix-vector product, the overall sensing
matrix $(\bm{\Lambda}~{}{\bf 0})+\tilde{\bm{E}}_{i}$ is dominated by the
diagonal matrix component. That is, $(\bm{\Lambda}~{}{\bf
0})+\tilde{\bm{E}}_{i}\approx(\bm{\Lambda}~{}{\bf 0})$. However, this implies
that the elements of ${\bm{a}}_{i}$ over the coordinates $[n]\setminus[K]$
cannot be recovered from ${\bm{y}}_{i}$ as the rows of the sensing matrix
supported only on $[K]$. It is thus necessary to exploit extra information to
recover ${\bm{a}}_{i}$.

In light of this, we note that $[{\bm{a}}_{i}]_{j}$ is zero for all $j$ in
${\mathcal{}S}_{i}$. The first matrix-vector product is thus an
$K-|{\mathcal{}S}_{i}|$-sparse vector which is supported on the set
${\mathcal{}S}_{i}^{c}=[K]-{\mathcal{}S}_{i}$. As we are interested in
studying a sufficient condition for perfect recovery, we see that (22) implies
the following linear system with the rows corresponding to
${\mathcal{}S}_{i}^{c}$ removed,

$\big{[}{\bm{y}}_{i}\big{]}_{{\mathcal{}S}_{i}}=\big{[}\tilde{\bm{E}}_{i}\big{]}_{{\mathcal{}S}_{i},:}{\bm{a}}_{i}\;.$
(23)

Compared to the original model (22), we can suppress the dominating diagonal
component in $\tilde{\bm{E}}_{i}$ using the support knowledge on
${\bm{a}}_{i}$.

The remaining task is to verify if the reduced sensing matrix
$\big{[}\tilde{\bm{E}}_{i}\big{]}_{{\mathcal{}S}_{i},:}$ is a good sensing
matrix. Notice that _dense_ and random matrices are known to exhibit good
properties in which one only requires $|{\mathcal{}S}_{i}|\geq
2\|{\bm{a}}_{i}\|_{0}\cdot\log n$ to achieve perfect recovery. On the other
hand, $\big{[}\tilde{\bm{E}}_{i}\big{]}_{{\mathcal{}S}_{i},:}$ is a _sparse_
matrix whose support depends on the out-neighbors of the nodes in
${\mathcal{}S}_{i}$. In particular, we have ${\rm
supp}(\big{[}\tilde{\bm{E}}_{i}\big{]}_{{\mathcal{}S}_{i},:})={\rm
supp}({\bm{A}}_{:,{\mathcal{}S}_{i}}^{\top})$. An interesting observation is
that the support of the sensing matrix depends on the support of ${\bm{A}}$
itself or the network that we wish to recover.

As it turns out, the perfect recovery condition for ${\bm{a}}_{i}$ boils down
to studying conditions on the _support_ of
$\big{[}\tilde{\bm{E}}_{i}\big{]}_{{\mathcal{}S}_{i},:}$. In particular, let
$G_{bi}(A,B)$ with $|A|=n$ and $|B|=|{\mathcal{}S}_{i}|\leq K$ be the bi-
partite graph representation of $[{\bm{E}}_{i}]_{{\mathcal{}S}_{i},:}$ such
that the adjacency matrix of $G_{bi}$, i.e.,
${\bm{A}}_{bi}\in\mathbb{R}^{|{\mathcal{}S}_{i}|\times n}$, has the same
support as $[{\bm{E}}_{i}]_{{\mathcal{}S}_{i},:}$.

###### Theorem 2

If $G_{bi}$ is an $(\alpha,\delta)$-unbalanced expander graph with left degree
bounded in $[d_{l},d_{u}]$ such that $2(d_{l}/d_{u})\cdot\delta>\sqrt{5}-1$
and $k\leq\frac{\alpha}{1+\rho\delta}n$, and $\tilde{\bm{E}}_{i}$ or
$-\tilde{\bm{E}}_{i}$ is non-negative, then the set

${\mathcal{}C}=\\{\hat{\bm{a}}_{i}~{}:~{}\big{[}\tilde{\bm{E}}_{i}\big{]}_{{\mathcal{}S}_{i},:}(\hat{\bm{a}}_{i}-{\bm{a}}_{i})={\bm{0}}\\}$
(24)

is a singleton for all $k$-sparse vector ${\bm{a}}_{i}$.

The proof of the theorem above can be found in subsection B.2. Consequently,
we observe that all the conditions in Theorem 2 hold, therefore solving (6)
yields $\hat{\bm{A}}={\bm{A}}$, i.e., we achieve perfect recovery of the
network.

###  B.1 A special case satisfying the conditions in Theorem 1

Let us first borrow the following proposition from [29]:

###### Proposition 2

[29, Proposition 5] Let $G(A,B)$ be a random bipartite with $|A|=n$,
$|B|=\beta\cdot n$ for some $\beta<1$ and the left degree bounded in
$[d_{l},d_{u}]$ (i.e., nodes in $A$ has bounded degree). Fix $\alpha\in(0,1)$
and the following holds:

$d_{l}\geq 3,~{}~{}\beta>\alpha,~{}~{}d_{l}>(H(\alpha)+\beta
H(\alpha/\beta))/(\alpha\log(\beta/\alpha))\;,$ (25)

where $H(\alpha)=\alpha\log\alpha+(1-\alpha)\log(1-\alpha)$ is the binary
entropy function. Then, with high probability as $n\rightarrow\infty$, the
graph $G(A,B)$ is an $(\alpha,1-1/d_{l})$-expander.

Following our discussions in Section 2.1, we suppose that the perturbed genes
$[K]$ are chosen such that for each $j\in V$, the $j$th gene has $d$ in-
neighbors in $[K]$. Notice that this requires selecting perturbation nodes
strategically, e.g., using advices from the experts.

Furthermore, the bipartite graph of the $K\times n$ matrix
$([{\bm{A}}]_{:,[K]})^{\top}$ may follow a regular random graph model with
constant degree $d$. Then, using similar arguments as in [29, Proposition 4],
it can be shown that the support of the sub-matrix
$([{\bm{A}}]_{:,{\mathcal{}S}_{i}})^{\top}$ corresponds to a random bipartite
graph with bounded degree in $[d-1,d]$ with high probability. The intuition is
that the matrix $([{\bm{A}}]_{:,{\mathcal{}S}_{i}})^{\top}$ is formed by
deleting a small number of rows from the $K\times n$ matrix
$([{\bm{A}}]_{:,[K]})^{\top}$.

Now fix $\alpha\in(0,1)$. Applying Proposition 2 shows that there exists
$\beta=c(d)\cdot\alpha$ with $c(d)>1$ such that this bipartite graph is an
$(\alpha,1-1/(d-1))$-expander with high probability and $\beta$ satisfying the
conditions in (25). Taking $|{\mathcal{}S}_{i}|=\beta n>c(d)\cdot\alpha n$, we
can substitute the numbers into the conditions (7) required by Theorem 1 to
obtain an upper bound on $\|{\bm{a}}_{i}\|_{0}$ where we can recover
${\bm{A}}$ perfectly. This yields the desirable result in (8).

###  B.2 Proof of Theorem 2

To prove Theorem 2, the following Lemma would be instrumental. Denote ${\rm
Null}({\bm{E}})$ as the null space of a matrix ${\bm{E}}$,
$S_{-}({\bm{w}})=\\{i\in[n]:w_{i}<0\\}$ and
$S_{+}({\bm{w}})=\\{i\in[n]:w_{i}>0\\}$ be the negative and positive support
of ${\bm{w}}$.

###### Lemma 1

If (i) $\overline{\bm{x}}$ is $k$-sparse, (ii) ${\bm{E}}\in\mathbb{R}^{m\times
n}$ satisfies that ${\bm{0}}\neq{\bm{w}}\in{\rm Null}({\bm{E}})$,
$|S_{-}({\bm{w}})|\geq k+1$, then the set
${\mathcal{}C}=\\{{\bm{x}}:{\bm{E}}({\bm{x}}-\overline{\bm{x}})={\bm{0}},~{}{\bm{x}}\geq{\bm{0}}\\}$
is a singleton.

Proof. Suppose $|{\mathcal{}C}|>1$ such that there exists
$\tilde{\bm{x}}\in{\mathcal{}C}$, $\tilde{\bm{x}}\neq\overline{\bm{x}}$. It is
straightforward to show that $\tilde{\bm{x}}=\overline{\bm{x}}+{\bm{w}}$,
where ${\bm{w}}\in{\rm Null}({\bm{E}})$. The assumption implies that
$|S_{-}({\bm{w}})|\geq k+1$ and
$\tilde{\bm{x}}=\overline{\bm{x}}+{\bm{w}}\ngeq{\bm{0}}$ as
$\overline{\bm{x}}$ is $k$-sparse. This contradicts
$\tilde{\bm{x}}\in{\mathcal{}C}$. Q.E.D.

The next step is to apply a generalization of [30, Theorem 4]:

###### Lemma 2

Let $n>m$ and ${\bm{E}}\in\mathbb{R}^{m\times n}$ be a non-negative matrix
that has the same support as the adjacency matrix of an
$(\alpha,\delta)$-unbalanced bipartite expander graph with left degrees
bounded in $[d_{l},d_{u}]$ and $\rho=d_{l}/d_{u}$. If
$\rho\delta>(\sqrt{5}-1)/2$, then for all $k$-sparse vector
$\overline{\bm{x}}$. The set

${\mathcal{}C}=\\{{\bm{x}}:{\bm{E}}({\bm{x}}-\overline{\bm{x}})={\bm{0}},~{}{\bm{x}}\geq{\bm{0}}\\},$
(26)

is a singleton if $k\leq\frac{\alpha}{1+\rho\delta}n$.

Proof. Using Lemma 1, it suffices to prove that any ${\bm{w}}\in{\rm
Null}({\bm{E}})$ with $S_{-}({\bm{w}})\subseteq A$, we have
$|S_{-}({\bm{w}})|\geq k+1$. We shall proceed by contradiction. Suppose that
there exists ${\bm{w}}\in{\rm Null}({\bm{E}})$ such that
$|S_{-}({\bm{w}})|\leq k$. Since $|S_{-}({\bm{w}})|\leq k\leq\alpha n$, the
expander property implies:

$\delta
d_{l}\cdot|S_{-}({\bm{w}})|\leq\delta|E(S_{-}({\bm{w}}),B)|\leq|N(S_{-}({\bm{w}}))|\;,$
(27)

where $N(S)$ denotes the set of neighbors of $S$ and $E(A,B)$ is the set of
edges between $A$ and $B$. The right hand side can be further upper bounded
as:

$|N(S_{-}({\bm{w}}))|\leq|E(S_{-}({\bm{w}}),B)|\leq
d_{u}\cdot|S_{-}({\bm{w}})|$ (28)

Moreover, we know that
$N(S_{-}({\bm{w}}))=N(S_{+}({\bm{w}}))=N(S_{-}({\bm{w}})\cup
S_{+}({\bm{w}}))$. This is because ${\bm{E}}{\bm{w}}={\bm{0}}$ and ${\bm{A}}$
is non-negative, thus the neighborhood sets must coincide to enforce nullity,
otherwise this will result in ${\bm{E}}{\bm{w}}\neq{\bf 0}$. Note that in
[30], the proof is achieved by assuming that ${\bm{E}}$ is the binary
adjacency matrix. We extend the same argument to the case when ${\bm{E}}$ is
non-negative. Using the above discussions and the inequality (27) applied on
$N(S_{+}({\bm{w}}))=N(S_{-}({\bm{w}}))$, we have:

$|S_{+}({\bm{w}})|\geq|N(S_{+}({\bm{w}})|/d_{u}\geq(d_{l}/d_{u})\delta|S_{-}({\bm{w}})|=\rho\delta|S_{-}({\bm{w}})|\;.$
(29)

As $S_{-}({\bm{w}})$ and $S_{+}({\bm{w}})$ are disjoint, we have
$|S_{+}({\bm{w}})\cup|S_{-}({\bm{w}})|\geq(1+\rho\delta)|S_{-}({\bm{w}})|$. As
such, we choose an arbitrary subset $\tilde{S}\subseteq S_{+}({\bm{w}})\cup
S_{-}({\bm{w}})$ such that $|\tilde{S}|=(1+\rho\delta)|S_{-}({\bm{w}})|$.
Notice that $|\tilde{S}|=(1+\rho\delta)|S_{-}({\bm{w}})|\leq\alpha n$. Using
the expander property again gives:

$|N(\tilde{S})|\geq\delta|E(\tilde{S},B)|\geq
d_{l}\delta(1+\rho\delta)|S_{-}({\bm{w}})|>d_{u}|S_{-}({\bm{w}})|$ (30)

As $\rho\delta>(\sqrt{5}-1)/2$, the last inequality is valid as
$\rho\delta(1+\rho\delta)>1$.

Finally, we reach a contradiction as

$|N(\tilde{S})|\leq|N(S_{-}({\bm{w}})\cup S_{+}({\bm{w}}))|=|N(S_{-}({\bf
w})|\leq d_{u}|S_{-}({\bm{w}})|,$ (31)

leading to $d_{u}|S_{-}({\bm{w}})|>d_{u}|S_{-}({\bm{w}})|$. The lemma is thus
proven. Q.E.D.

Identifying that $G_{bi}$ satisfies the conditions in Lemma 2, our desirable
results can be obtained.

##  Appendix C Derivation of Problem (12)

From (11), for each $i\in[n]$, we can model the vector ${\bm{y}}_{i}$ as

${\bm{y}}_{i}={\bm{H}}_{i}{\bm{a}}_{i}+\bm{\Delta}_{i}{\bm{a}}_{i}+\bm{\epsilon}\;,$
(32)

where $\bm{\Delta}_{i}$ models the error in the observed measurement matrix
${\bm{H}}_{i}$. Let $r>0$, the uncertainty set for $\bm{\Delta}_{i}$ is
defined such that each row vector in the matrix has a bounded norm, i.e.,

${\mathcal{}U}_{r}=\\{\bm{\Delta}_{i}~{}:~{}\|{\bm{d}}_{k}\|_{2}\leq
r,~{}\forall~{}k\in[K+1]\\}\;,$ (33)

where ${\bm{d}}_{k}$ is the $k$th row vector of $\bm{\Delta}_{i}$.

To recover ${\bm{a}}_{i}$, we consider minimizing the following robust
objective:

$J(\hat{\bm{a}}_{i})=\lambda\|\hat{\bm{a}}_{i}\|_{1}+\max_{\bm{\Delta}_{i}\in{\mathcal{}U}_{r}}\|{\bm{y}}_{i}-{\bm{H}}_{i}\hat{\bm{a}}_{i}-\bm{\Delta}_{i}\hat{\bm{a}}_{i}\|_{2}\;,$
(34)

which can be upper bounded by:

$\begin{split}J(\hat{\bm{a}}_{i})&\leq\lambda\|\hat{\bm{a}}_{i}\|_{1}+\|{\bm{y}}_{i}-{\bm{H}}_{i}\hat{\bm{a}}_{i}\|_{2}+\max_{\bm{\Delta}_{i}\in{\mathcal{}U}_{r}}\|\bm{\Delta}_{i}\hat{\bm{a}}_{i}\|_{2}\\\
&=\lambda\|\hat{\bm{a}}_{i}\|_{1}+\|{\bm{y}}_{i}-{\bm{H}}_{i}\hat{\bm{a}}_{i}\|_{2}+r\sqrt{K+1}\cdot\|\hat{\bm{a}}_{i}\|_{2}\;,\end{split}$
(35)

where the last equality is achieved by applying Cauchy-Schwarz and setting
each row of $\bm{\Delta}_{i}$ to
$r\cdot\hat{\bm{a}}_{i}/\|\hat{\bm{a}}_{i}\|_{2}$. Setting
$\gamma=r\sqrt{K+1}$ and minimizing the upper bound function yields the robust
network recovery problem (12).

##  Appendix D Details of experiments on empirical data

This section describes the pre/post-processing procedures for the empirical
support data.

Subspace projection for ‘denoising’. We apply a subspace projection method as
a pre-processing stage in the RIDS method. In particular, let the set of chips
(or vectors of expression levels) taken under the _no perturbation_ condition
be ${\mathcal{}C}_{nopert}$, we can write the gene expression levels from the
$c$th chip as

$\overline{\bm{x}}_{obs,c}=\overline{\bm{x}}+{\bm{w}}_{c},~{}\forall~{}c\in{\mathcal{}C}_{nopert}\Longrightarrow\underbrace{\big{(}\ldots~{}\overline{\bm{x}}_{obs,c}~{}\ldots\big{)}}_{\mathrel{\mathop{:}}=\overline{\bm{X}}_{obs}}=\overline{\bm{x}}{\bf
1}^{\top}+{\bm{W}}_{obs}\;,$ (36)

where $\overline{\bm{x}}$ is the unperturbed steady state satisfying
$\overline{\bm{x}}={\bm{A}}{\bm{h}}(\overline{\bm{x}})$ (cf. Eq. (2)) and
${\bm{w}}_{c}$ is modeled as an additive noise. When the noise is small, we
observe that the $n\times|{\mathcal{}C}_{nopert}|$ matrix
$\overline{\bm{X}}_{obs}$ formed by stacking up $\overline{\bm{x}}_{obs,c}$
horizontally is close to rank-one. In light of this, we recover
$\overline{\bm{x}}$ by taking the top left-singular vector of
$\overline{\bm{X}}_{obs}$, i.e.,

$\tilde{\overline{\bm{x}}}=\sigma_{1}(\overline{\bm{X}}_{obs})\cdot{\bm{u}}_{1}~{}~{}~{}\text{where}~{}~{}~{}\overline{\bm{X}}_{obs}={\bm{U}}\bm{\Sigma}{\bm{V}}^{\top}\;,$
(37)

and $\sigma_{1}(\overline{\bm{X}}_{obs})$ is the largest singular value of
$\overline{\bm{X}}_{obs}$. Notice that $\hat{\overline{\bm{x}}}$ lies in the
null space of the sub-matrix $[{\bm{U}}]_{:,\text{2:end}}$. For each
$k\in[K]$, we can estimate the difference ${\bm{x}}[k]-\overline{\bm{x}}$ by

$\widetilde{(\bm{x}[k]-\overline{\bm{x}})}=\big{(}[{\bm{U}}]_{:,\text{2:end}}[{\bm{U}}]_{:,\text{2:end}}^{\top}\big{)}\cdot\frac{1}{|{\mathcal{}C}_{per}[k]|}\sum_{c\in{\mathcal{}C}_{per}[k]}\bm{x}_{obs,c}\;,$
(38)

where $\bm{x}_{obs,c}$ is the gene expression levels from chip $c$ and
${\mathcal{}C}_{per}[k]$ is the set of chips that corresponds to the $k$th
perturbation experiment. Naturally, we have
$\tilde{\bm{x}}[k]=(\widetilde{\bm{x}[k]-\overline{\bm{x}}})+\tilde{\overline{\bm{x}}}$.

Normalizing the gene expression data. The vectors of gene expression level
$\tilde{\overline{\bm{x}}}$ and $\tilde{\bm{x}}[k]$ obtained from from the
above are normalized through dividing the vectors by the constant

$c_{norm}\mathrel{\mathop{:}}=\frac{1}{2}\max\Big{\\{}\max_{i\in[n]}~{}(\tilde{\overline{x}}_{i}),~{}\max_{k\in[K]}\max_{i\in[n]}~{}(\tilde{x}_{i}[k])\Big{\\}}\;$
(39)

such that the values of the normalized gene expression level ranges in
$[0,2]$.

Parameters of RIDS. For the zeros index set ${\mathcal{}S}$ in (5), we set
$\delta=0.005$. The regularization parameters $\rho$ and $\gamma$ are set to
$3\times 10^{0}$ and $5\times 10^{0}$ (resp. $3\times 10^{-1}$ and $5\times
10^{-1}$) in (12) for _E. coli_ (resp.  _S. cerevisiae_) network. Meanwhile,
the AO procedure (13) aims at tackling (12) with the model response function
template (14), which has several parameters to be initialized — for each
$i\in[n]$, we initialize the AO procedure by setting $b_{1}=0.5$ and
$b_{2}=0.5$ and the step size $\epsilon$ is set to $1\times 10^{-2}$. Lastly,
the AO procedure terminates after $L=10$ iterations in the numerical
experiment for a better trade-off between complexity and accuracy.

Post-processing of the estimated GRN. The proposed algorithm may recover also
the negative edge weights in $\hat{\bm{A}}$. However, we shall treat these
negative edge weights as positive ones for the ease of benchmarking using the
DREAM5 evaluation script. In particular, we consider $|\hat{\bm{A}}|$ as our
estimate of the GRN. Furthermore, we normalize the elements in
$|\hat{\bm{A}}|$ by $\max_{i,j}|[\hat{\bm{A}}]_{ij}|$ such that they range in
$[0,1]$. For the AUROC/AUPR evaluation, we only count the top 100,000 (and
500,000) predicted links made by our algorithm.

##  Appendix E Additional Results on Empirical Data

We present additional experiment results for the _in vivo_ dataset from
DREAM5.

Figure 5: ROC and Precision-recall curves for the _in vivo_ networks. (Top)
_E. coli_ network. (Bottom) _S. cerevisiae_ network. These curves are produced
using the DREAM5 evaluation script and the AUROC/AUPR values can be found in
Table 1.

###  E.1 GRN recovery without the aid of transcription factors

For simplicity, we shall only focus on the _E. coli_ network. We consider
applying the RIDS method without the help of the information provided from the
list of transcription factors. In particular, we solve (12) by fixing the
model parameter at ${\bm{b}}=(0.047,0.5893)$ and set $\rho=3$, $\gamma=5$.

AUROC and AUPR scores. Even without any aid from the list of transcription
factors, the RIDS method continues to yield a top AUROC/AUPR performance. In
particular, the top $100,000$ predicted links from the identified
$\hat{\bm{A}}$, restricted to the DREAM5 identified transcription factors,
yields an AUROC of $0.6745$ and AUPR of $0.0540$ for the _E. coli_ network, as
seen in Table 1.

Figure 6: Inferring directionality of the GRN using RIDS — comparison of the
inferred regulatory strengths and the received regulatory strengths for each
gene in _E. coli_. The top figure shows the regulatory strengths for each gene
exerted on the other genes; the bottom figure shows the regulatory strengths
_received_ for each gene from the other genes.

Directionality of the GRN. We test whether the directionality in GRNs can be
inferred under the scenario described in this subsection i.e., without using
the provided list of transcription factors (TFs). To emphasize on the links
with higher weights, we normalize the inferred GRN $\hat{\bm{A}}$ as
$[\widetilde{\bm{A}}]_{ij}\mathrel{\mathop{:}}=|\hat{A}_{ij}|^{q}/\sum_{i^{\prime},j^{\prime}}|\hat{A}_{i^{\prime}j^{\prime}}|^{q}$
for $q=2$. Then we compare the following for each gene $i\in[n]$ —

  * •

Regulatory strengths — absolute sum of the columns of $\widetilde{\bm{A}}$,
i.e., $\sum_{j=1}^{n}|\widetilde{A}_{ji}|$.

  * •

Strengths to be regulated — absolute sum of the rows of $\widetilde{\bm{A}}$,
i.e., $\sum_{j=1}^{n}|\widetilde{A}_{ij}|$.

The comparison can be found in Figure 6. We observe that the inferred GRN is
directional as fewer genes have significantly higher regulatory strengths
exerted on the others, while the regulatory strengths received for the genes
are more uniformly distributed.