\section*{Introduction}
Complex interactions between many genes give rise to the biological
structure and function that sustain life. The Central Dogma [\citet
{jacob61,crick70}]\vadjust{\goodbreak} provides a qualitative description of how these
processes occur, but precise quantitative modeling is still needed
[\citet{tyson03,rosenfeld11}]. Research into the detailed mechanisms of
gene expression over the past few decades has shown that expression is
regulated by a complex system of gene interactions. Recently,
microarray and sequencing technologies [\citet{derisi97,ren00,robertson07,mortazavi08}] have enabled high-throughput genome-wide
expression level measurements. This data enables detailed study of gene
networks [\citet{holstege98,lee02,tegner03,segal03,bar-joseph03,hu07,zhou07}]. The goal is to understand how genes interact to give
rise to the biochemical complexity that allows organisms to live, grow
and reproduce.

Gene expression measurements contain information useful for
reconstructing the underlying interaction structure [\citet{derisi97,holstege98,hughes00}] because gene regulatory systems have a defined
ordering [\citet{avery92}], forming pathways that connect to form
networks [\citet{alon07,desmet10}]. Many gene regulation pathways have
been discovered over the past few decades [\citet{hartwell10,alberts07}]. At the turn of the century, researchers began applying
statistical tools to genome-wide expression data to understand complex
gene interactions. Eisen et al. showed that genes from the same
pathways and with similar functions cluster together by expression
pattern [\citet{eisen98}]. Soon afterward, module-based network inference methods
appeared, which group co-expressed genes into cellular function modules
[\citet{segal03,bar-joseph03}]. Recently, methods based on descriptive
but nonmechanistic mathematical models [\citet{gardner03,tegner03,bansal07,faith07,friedman04}] have gained prominence. These models
describe gene regulation quantitatively and can be used to simulate and
predict systems behaviors [\citet{palsson11,dehmer11}]. However, more
work is needed to develop effective model-based methods for inferring
gene network structure from experimental data.

Existing inference methods typically rely either on heuristic
approaches or on very simple, local models, like linear differential
equation models in a neighborhood of a particular steady-state.
Statistical corrrelation is a common method of establishing network
connections [\citet{dehmer11}] and can be very useful when hundreds or
thousands of genes are monitored under specific, local cellular
conditions (e.g., for grouping genes with similar functions).
However, this approach works poorly when perturbations drive the
network far from the original steady-state. Global nonlinear models are
essential to account for complex global system behaviors, like the
transformation of a normal cell into a cancerous cell due to the
amplification of a particular gene.

As a basis for our inference approach, we chose a standard global
nonlinear model: the quantitative, experimentally interpretable
biophysics-based ordinary differential equation (ODE) gene regulation
model of Bintu et al. (\citeyear{bintu05models,bintu05apps}). Many
models of this type have been proposed, and the idea traces back to the
beginning of systems biology in the biophysics field
[\citet{ackers82,shea85,vonhippel74}], but the Bintu model is
widely accepted within the biophysics community
[\citet{bintu05models}]. The Bintu model is based on the
thermodynamics of RNA transcription, the process at the core of gene
expression regulation [\citet{holstege98,hu07}]. Transcription
occurs when RNA polymerase (RNAP) binds the gene promoter;
transcription factors (TFs) can modulate the RNAP binding energy to
activate or repress transcription. RNA transcripts are then translated
into protein. Bintu models the mechanism of transcription in detail,
using physically interpretable parameters. The form of the equations is
rich and flexible enough to include the full range of gene regulatory
behavior. Another notable biophysics-based model is that of the annual
DREAM competition, but it has many biochemical assumptions and model
parameters, like the Hill coefficient of transcription factor binding
events, that cannot be estimated using gene expression measurements, so
the network reconstruction requires ad hoc inference methods to learn
the underlying gene interactions
[\citet{yip10,pinna10,marbach10,schaffter11}]. Compared to the
DREAM model, the Bintu model has the advantages of simplicity and
interpretability, and better lends itself to principled inference.

In this paper, we propose an experimental design and associated
statistical method for inferring an unknown gene network by fitting the
ODE-based Bintu gene regulation model. The required data is gene
expression measurements at a set of perturbed steady-states induced by
gene knockdown and overexpression [\citet{Huang05}]. We show how to
design a sequence of experiments to collect the data and how to use it
to fit the parameters of the Bintu model, leading to a set of ODEs that
quantitatively characterize the regulatory network. Although the
original fitting problem is nonlinear, we can transform it into a
convex optimization problem by restricting our attention to
steady-states. We use the lasso [\citet{tibs96}] for parameter
selection. As a proof of principle, we test the method on a simulated
embryonic stem cell (ESC) transcription network [\citet{chickarmane08}]
given by a system of ODEs based on the Bintu model. Here, we
demonstrate that the inference algorithm is computationally efficient,
accounts for synergistic regulation and self-regulation, and correctly
recovers the parameters used to generate the data. Furthermore, the
method requires only a set of steady-state gene expression
measurements. Experimental researchers in the biological sciences can
use this method to infer gene networks in a much more principled,
detailed manner than earlier approaches allowed.

\section*{Dynamical systems model}
We model gene expression regulation as a dynamical system. Let $x \in\R
^n$ represent RNA concentrations and $y \in\R^n$ represent protein
concentrations corresponding to a set of $n$ genes. We assume that the
production rate of the RNA transcript $x_i$ of gene $i$ is proportional
to the probability $f(y)$ that RNA polymerase (RNAP) is bound to the
promoter. That is, we assume that RNA transcription occurs at a rate
$\tau_i$ whenever RNAP is bound to the promoter. We model the
probability that RNAP is bound to the promoter as a nonlinear function
$f$ of $y$, since RNAP binding is regulated by a set of TFs. Further,
we assume that the production rate of protein product $y_i$ of gene $i$
is proportional to the concentration of the RNA transcript $x_i$, and
that both the RNA transcript and protein products of gene $i$ degrade
at fixed rates ($\lambda_i^{\mathrm{RNA}}$, $\lambda_i^{\mathrm{Protein}}$),
\begin{eqnarray}\label{eqexpression}
\frac{dx_i}{dt} &=& \tau_i f_i(y) -
\lambda_i^{\mathrm{RNA}} x_i,
\nonumber\\[-8pt]\\[-8pt]
\frac{dy_i}{dt} &=& r_ix_i - \lambda_i^{\mathrm{Protein}}
y_i.\nonumber
\end{eqnarray}
Based on the thermodynamics of RNAP and TF binding, one can deduce the
following form for $f_i$ [Bintu et~al. (\citeyear{bintu05models,bintu05apps})]:
\begin{equation} \label{eqf}
f_i(y) = \frac{b_{i0} + \sum_{j=1}^m b_{ij} \Pi_{k \in S_{ij}} y_k}{1
+ \sum_{j=1}^m c_{ij} \Pi_{k \in S_{ij}} y_k},
\end{equation}
where $S_{ij}$ lists the gene products that interact to form a
regulatory complex, and $b_{ij}, c_{ij}$ are nonnegative coefficients
that must satisfy $c_{ij} \ge b_{ij} \ge0$. (We assume that the
concentration of each complex is proportional to the product of the
concentrations of the constituent proteins, and absorb the
proportionality constant into corresponding coefficients $b_{ij},
c_{ij}$.) The coefficients $b_{ij}$ and $c_{ij}$ depend on the binding
energies of regulator complexes to the promoter. $b_{i0}$ and $c_{i0}$
correspond to the case when the promoter is not bound by any regulator
($\Pi_{k \in S_{i0}} y_k = 1$), and the coefficients are normalized so
that $c_{i0} = 1$. Details and a derivation are given in the \hyperref[app]{Appendix}.

The form of $f_i$ allows us to model the full spectrum of regulatory
behavior in quantitative detail. Terms that appear in the denominator
only are repressors, and the degree of repression depends on the
magnitude of the coefficient, while terms that appear in the numerator
and denominator may act as either activators or repressors depending on
the relative magnitudes of the coefficients and the current gene
expression levels. Terms may represent either single genes or gene
complexes. The model can even be extended to account for environmental
factors that affect gene regulation, though we will not discuss it
further here.

\begin{figure}

\includegraphics{645f01.eps}

\caption{Simple two-gene network example described by equation
(\protect\ref{eqmodelexample}) (with parameters $b_{11} = c_{11} =
0.1$ for activators; $c_{12} = 10$ for repressors; and $b_{10} = 0.01$
for constants in the numerator). Gene 1 is activated by the protein
product of gene 2 and repressed by its own product (an example of
self-regulation). Gene 2 is repressed by a complex formed by the
product of gene 1 and its own product (synergistic self-regulation). In
the diagram, the edge colors indicate activation (green) or repression
(red) and the edge weights indicate coefficient sizes, illustrated
above with typical sizes.} \label{figmodelexample}
\end{figure}

As an example, consider the simple two-gene network shown in
Figure~\ref{figmodelexample}. Suppose that genes 1 and 2 have RNA
concentrations $x_1, x_2$, and protein concentrations $y_1, y_2$,
respectively, and that gene 1 is activated by protein 2 and repressed
by its own product (protein 1), while gene 2 is repressed by a complex
formed by proteins 1 and 2. The situation corresponds to the following
equations:
\begin{eqnarray} \label{eqmodelexample}
\frac{dx_1}{dt} &=& \tau_1 \frac{b_{10} + b_{11} y_2 }{1 + c_{11}
y_2 + c_{12} y_1} -
\lambda_1^{\mathrm{RNA}} x_1,\qquad \frac{dy_1}{dt} =
r_1x_1 - \lambda_1^{\mathrm{Protein}}
y_1,
\nonumber\\[-8pt]\\[-8pt]
\frac{dx_2}{dt} &=& \tau_2 \frac{b_{20} }{1 + c_{21} y_1 y_2} - \lambda
_2^{\mathrm{RNA}} x_2,\qquad \frac{dy_2}{dt} =
r_2x_2 - \lambda_2^{\mathrm{Protein}}
y_2.\nonumber
\end{eqnarray}
In the notation above, we have $S_{11} = \{ 2 \}, S_{12} = \{ 1 \},
S_{21} = \{ 1, 2 \}$. The parameters $b_{10}, b_{11}, c_{11},\ldots$
determine the magnitude of the repression or activation. As this
example shows, the model is flexible enough to capture a wide range of
effects, including self-regulation (i.e., regulation of a gene by
its own protein product, most commonly as repression) and synergistic
regulation by protein complexes (two or more proteins bound together to
form a regulatory unit), in quantitative detail. Furthermore, the model
is predictive: if we know or can infer the coefficients in the model,
we can predict the future behavior of the system starting from any
initial condition.

\section*{Inference problem}

The model given by equations (\ref{eqexpression}) and (\ref{eqf}) fully
describes the evolution of RNA and protein levels and provides a
comprehensive, quantitative model of gene regulation, provided we know
the parameters. Unfortunately, $b_{ij},c_{ij}$ are extremely difficult
to measure, as they depend on binding energies of RNAP and TFs to the
gene promoter. The sheer number of measurements required to
characterize all possible TFs (both individual proteins and complexes)
also makes this approach infeasible. Therefore, our goal is to use a
systems level approach to fit the model using RNA expression data.
Specifically, we will assume that $\tau_i, \lambda_i^{\mathrm{RNA}}, \lambda
_i^{\mathrm{Protein}}$ are known or can be measured (if these quantities are not
available, we can simply absorb them into the coefficients
$b_{ij},c_{ij}$, although more accurate rate estimates will likely
improve the coefficient estimates). Our data will be measurements of
the RNA concentrations $x$ at many different cellular steady-states
(which correspond to steady-states of the dynamical system). The
problem is to infer the values of the coefficients $b_{ij},c_{ij}$.

\section*{Linear problem at steady-state}

The key to solving this problem efficiently is to restrict our
attention to steady-states, as proposed by Choi (\citeyear{choi12}). This
restriction allows us to transform a nonlinear ODE fitting problem into
a linear regression problem. A steady-state of the system is one in
which RNA and protein levels are constant: $\frac{dx_i}{dt} = \frac
{dy_i}{dt} = 0$. Steady-states of the system correspond to cell states
with roughly constant gene expression levels, like embyronic stem cell,
skin cell or liver cell. In contrast, an embryonic stem cell in the
process of differentiating is not in steady-state. Perturbed
steady-states are particularly interesting. After a perturbation like
gene knockdown, a cell's gene expression levels are in flux for some
time while they adjust to the change. Eventually, if it is still
viable, the cell may settle to a new steady-state [\citet{Huang05}].
These perturbed steady-states are especially helpful for understanding
gene regulation.

In our model, the steady-state conditions $\frac{dx_i}{dt} = \frac
{dy_i}{dt} = 0$ mean that
\[
0 = \tau_i f_i(y) - \lambda_i^{\mathrm{RNA}}
x_i,\qquad 0 = r_ix_i - \lambda_i^{\mathrm{Protein}}
y_i \quad\implies\quad y_i = \frac{r_i x_i}{\lambda
_i^{\mathrm{Protein}}}.
\]
Defining $\tilde{f}_i(z) = f_i(\frac{r_i}{\lambda_i^{\mathrm{Protein}}} z) $ yields
\[
0 = \tau_i \tilde{f}_i(x) - \lambda_i^{\mathrm{RNA}}x_i.
\]
Absorbing the constants into the coefficients $b_{ij},c_{ij}$ (so that
$\tilde{b}_{ij} =\break b_{ij} \Pi_{k \in S_{ij}} \frac{r_k}{\lambda
_k^{\mathrm{Protein}}}$, $\tilde{c}_{ij} = c_{ij} \Pi_{k \in S_{ij}} \frac
{r_k}{\lambda_k^{\mathrm{Protein}}}$), we obtain the final equation
\[
\tau_i \frac{b_{i0} + \sum_j b_{ij} \Pi_{k \in S_{ij}} x_k}{1 + \sum_j
c_{ij} \Pi_{k \in S_{ij}} x_k} - \gamma_i x_i
= 0
\]
or
\[
\tau_i \biggl(b_{i0} + \sum_j
b_{ij} \Pi_{k \in S_{ij}} x_k\biggr) -
\gamma_i x_i \biggl(1 + \sum
_j c_{ij} \Pi_{k \in S_{ij}} x_k
\biggr) = 0
\]
(by multiplying both sides by the denominator). The last equation is
linear in the coefficients $b_{ij},c_{ij}$!
In order to solve for $b_{ij},c_{ij}$, we will need to collect many
different expression measurements $x$ at both naturally occurring and
perturbed steady-states. Each steady-state measurement will lead to a
different linear equation. These equations can be arranged into a
linear system that we can solve for the coefficients.

\section*{Problem formulation}
Our problem is to find $b_{ij},c_{ij}$ such that
\begin{eqnarray*}
0 &=& \tau_i \biggl(b_{i0} + \sum
_j b_{ij} \Pi_{k \in S_{ij}} x_k^{(m)}
\biggr) \\
&&{}- \gamma_i x_i \biggl(1 + \sum
_j c_{ij} \Pi_{k \in S_{ij}} x_k^{(m)}
\biggr)\qquad \forall m = 1,\ldots, M,
\end{eqnarray*}
given RNA expression data $x^{(m)}$ at many different steady-state
points $m=1,\ldots,M$ and known translation and degradation rates $\tau
_i, \lambda_i^{\mathrm{RNA}}, \lambda_i^{\mathrm{Protein}}$. (The experimental means of
collecting the necessary steady-state expression data will be discussed
in the next section.) We solve a separate problem for each gene $i$,
since the coefficients $b_{ij}, c_{ij}$ in the differential equation
$dx_i/dt = \cdots$ for gene $i$ are independent of the coefficients in
the differential equations for other genes. Since we cannot know ahead
of time which potential regulatory terms $\Pi_{k \in S_{ij}} x_k$ are
actually involved, we include all possible terms up to second-order and
look for sparse $b_{ij},c_{ij}$, intepreting $c_{ij} = 0$ to mean that
term $\Pi_{k \in S_{ij}} x_k$ is not a regulator of gene $i$.

Consider gene 2 in the two-gene example. Suppose we have expression
measurements for a naturally occurring steady-state $(x_1^0, x_2^0)$,
and a perturbed steady-state following gene 1-knockout $(0, x_2^1)$. We
obtain two linear equations in the coefficients $b_{20}, c_{21}$:
\begin{eqnarray*}
\tau b_{20} - \lambda_2 x_2^0
\bigl(1 + c_{21} x_1^0 x_2^0
\bigr) &=& 0\qquad \bigl(\mbox{steady-state } \bigl(x_1^0,
x_2^0\bigr)\bigr),
\\
\tau b_{20} - \lambda_2 x_2^1 &=&
0\qquad \bigl(\mbox{steady-state } \bigl(0, x_2^1\bigr)\bigr).
\end{eqnarray*}
If we knew a priori that complex $x_1 x_2$ was the only regulator of
gene 2, these two equations would allow us to solve for the
coefficients ($b_{20} = \frac{\lambda_2 x_2^1}{\tau}$, $c_{21} = \frac
{x_2^0 - x_2^1}{(x_2^0)^2}$). Typically we do not know the regulators
beforehand, however, and we need to use the data to identify them. That
is, we include all possible terms (up to second-order) in the equations:
\begin{eqnarray*}
&&
\tau\bigl(b_{20} + b_{21}x_1^{(m)}
x_2^{(m)} + b_{22} x_1^{(m)}
+ b_{23} x_2^{(m)}\bigr)\\
&&\quad{} - \lambda_2
x_2^{(m)} \bigl(1 + c_{21}x_1^{(m)}
x_2^{(m)} + c_{22} x_1^{(m)}+
c_{22} x_2^{(m)}\bigr) \\
&&\qquad= 0
\end{eqnarray*}
and estimate sparse coefficients $b_{ij}, c_{ij}$ using several
steady-state measurements $(x_1^{(m)}, x_2^{(m)})$. (We should find
that the recovered coefficients $b_{21}, b_{22},\break b_{23}, c_{22},
c_{23}$ are very close to zero, since the corresponding terms do not
appear in the true equation.)

Temporarily suppressing the superscript $m$ denoting the observation,
we can compactly express the general system above by defining $z_i$ as
the vector with entries $z_{i}(j) = \Pi_{k \in S_{ij}} x_k$ [with the\vadjust{\goodbreak}
convention\vspace*{1pt} that $z_i(0) = 1$, $z_i(j) = x_j$ for $j=1,\ldots,n$], which yields
\[
0 = \tau_i b_i^Tz_i -
\gamma_i z_i(i) c_i^Tz_i
\]
for each observation ($m=1,\ldots, M$). If we form a matrix $G_i$ by
concatenating the row vectors $z_i^{(1)},\ldots, z_i^{(M)}$ and let
$D_i$ be a diagonal matrix with entries $z_i^{(m)}(i),   m =
1,\ldots,M$, we can express this as
\[
\left[\matrix{\tau_i G_i & -
\gamma_i D_i G_i}\right] \left[\matrix{b_i
\cr
c_i}\right] = 0
\]
with the constraints $0 \le b_i \le c_i,   c_i(0) = 1$. Stating the
problem in this form elucidates the required number of steady-state
measurements, $M$. If the linear system above were dense and had no
constraints on the coefficients $b_{ij}, c_{ij}$, and the steady-state
expression vectors were (numerically) linearly independent, then we
would require $M = 2 T_{n,k}$, where $T_{n,k}$ is the number of the
terms in the rational-form polynomial of degree $k$ in $n$ genes ($k=2$
if we include up to second-order regulatory interactions). $T_{n,k}$ is
equal to the number of subsets of $\{1,2,\ldots,n\}$ with $k$ or fewer
elements since each term represents an interaction between $j$ distinct
genes ($0 \le j \le k$), hence, $T_{n,k} = \sum_{j=0}^k { n \choose j }
\le n^k $ for $k \le n$ ($T_{n,2} \le n^2$, e.g.). However, the
constraints reduce the dimension of the solution space [$c_i(0) = 1$
reduces it by 1, while $0 \le b_i \le c_i$ reduces it by up to $n$],
and our algorithm also uses $\ell_1$-regression to search for sparse
solutions, which may allow us to reconstruct the coeffcients from far
fewer measurements than $2 T_{n,k}$.

\section*{Experimental approach}

The set of steady-state gene expression measurements needed to fit the
model can be generated via a systematic sequence of gene perturbation
experiments. Figure~\ref{figexperiments} summarizes the overall
approach to finding the regulatory interactions among a set of genes
\begin{figure}

\includegraphics{645f02.eps}

\caption{Experimental approach for gene network inference. (1) Design
and perform perturbation experiments targeting each gene (or possibly
pair of genes) in the network: these may include overexpression,
knockdowns or knockouts. (2) Following each perturbation, allow the
system to settle to a new steady-state. (3) Measure expression levels
of all genes at each induced steady-state, and collect results in a
data matrix. (4) Use steady-state expression data as input to inference
algorithm. (5) Construct regulatory network from inference algorithm output.}
\label{figexperiments}
\end{figure}
comprising a (roughly) self-contained network of interest. First,
molecular perturbations targeting each gene, or possibly pair of genes,
in the network would be designed and applied one at a time. Following
each perturbation, the cells would be allowed to settle down to a new
steady-state, at which point the gene expression levels would be
measured. The collection of gene expression measurements from different
steady-states would be input to the inference algorithm described in
the next section, which outputs a dynamical systems model of the gene
network capable of predicting the behavior of the network following
other perturbations. Perturbation data not used in the inference
algorithm could be used to validate the recovered model.

The key experimental steps in this procedure, gene perturbations and
gene expression measurements, are established technologies. Gene
perturbations, including overexpression, knockdown and knockout, are
routinely used in biological studies to\vadjust{\goodbreak} investigate gene function.
These experiments can be performed for many laboratory organisms and
cell lines both in vitro and in vivo [\citet{alberts07}].
Overexpression experiments amplify a gene's expression level, usually
by introducing an extra copy of the gene. Knockdown experiments
typically use RNAi technology: the cell is transfected with a short DNA
sequence, driven by a (possibly inducible) promoter element, that
produces siRNA or shRNA that specifically binds the RNA transcripts of
the gene of interest and triggers degradation. Morpholinos can also be
used for gene knockdown. Gene knockout can be achieved by removing all
or part of a gene to permanently disrupt transcription [\citet
{alberts07}]. Overexpression [\citet{rodriguez07}], knockdown [\citet
{rodriguez07,foygel08}] and knockout [\citet{lengner07}] experiments
have all been performed for the Oct4 gene, which helps maintain the
stem cell steady-state. In some cases,\vadjust{\goodbreak} much of the work is already
done: for example, the \emph{Saccharomyces} Genome Deletion Project has
a nearly complete library of deletion mutants [\citet{winzeler99}].

Techniques for gene expression measurement are also well-established.
Gene expression is usually measured at the transcript level: the RNA
transcripts are extracted and reverse-transcribed into cDNA, which can
be quantified with either RT-qPCR, microarray or sequencing
technologies [\citet{alberts07,mortazavi08}]. Housekeeping gene
expression measurements are used as controls to determine the
expression levels of the genes of interest. The gene perturbations and
subsequent expression measurements required to collect data for our
inference algorithm may be time-consuming due to the large number of
perturbations, but all the experimental techniques are quite standard
and resources like deletion libraries can be extremely helpful.

\section*{Algorithm}
We need to solve the linear system
\[
\left[\matrix{\tau_i G_i & -
\gamma_i D_i G_i}\right] \left[\matrix{b_i
\cr
c_i}\right] = 0
\]
for $b_i, c_i$, subject to the constraints $0 \le b_i \le c_i$, $c_i(0)
= 1$. To account for measurement noise and encourage sparsity in
$b_{i}, c_{i}$ (since we know that each gene has only a few
regulators), we will minimize the $\ell_2$-norm error with $\ell_1$
regularization [\citet{tibs96}], which leads to the convex optimization problem
\begin{eqnarray}\label{eqopt}
&&\mbox{minimize } \biggl\| \left[\matrix{\tau_i G_i & -
\gamma_i D_i G_i}\right] \left[\matrix{b_i
\cr
c_i}\right] \biggr\|_2 + \lambda\bigl(\| b_i
\|_1 + \| c_i \|_1\bigr)
\nonumber\\[-8pt]\\[-8pt]
&&\quad\mbox{subject to}\quad 0 \le b_i \le c_i,\qquad
c_i(0) = 1,\nonumber
\end{eqnarray}
where $\lambda$ is a parameter controlling sparsity that we can choose
using cross-validation. Since the problem is convex, it can be solved
very efficiently even for large values of $n$ and $m$.

\section*{Nonidentifiability}

Our model's ability to capture self-regulation is very powerful, but it
also leads to a particular form of nonidentifiability. For certain
forms of the equation, given only steady-state measurements, it can be
impossible to determine whether self-regulation is either completely
absent or present in every term. Specifically, any valid equation of
the form
\begin{equation}
\label{eqsimple} \frac{dx_i}{dt} = \frac{b_{i0} + \sum_{j=1}^N b_{ij} \Pi_{k \in
S_{ij}} x_k}{1 + \sum_{j=1}^N c_{ij} \Pi_{k \in S_{ij}} x_k} - \gamma_i
x_i,\qquad b_{i0} < 1,
\end{equation}
is indistinguishable at steady-state from any member of the following
family of valid equations indexed by the constant $w$:
\begin{equation}
\label{eqhighorder}\quad \frac{dx_i}{dt} = \frac{ (wb_{i0} + \gamma_i)x_i + \sum_{j=1}^N
wb_{ij} \Pi_{k \in{S}_{ij}} x_ix_k}{1 + wx_i + \sum_{j=1}^N wc_{ij} \Pi
_{k \in{S}_{ij}} x_ix_k} -
\gamma_i x_i,\qquad w \ge\frac{\gamma}{1 - b_{i0}}.
\end{equation}
We will refer to these as the ``simple'' and ``higher-order'' forms of the
equation, respectively.
The short proof of their equivalence is given in section S1
of the supplementary article [\citet{meister13supplement}]. The
condition $w \ge\frac{\gamma}{1 - b_{i0}}$ guarantees that $w > 0$ and
$0 \le wb_{i0} + \gamma_i \le w$ (since $0 \le b_{i0} < 1$) and $0 \le
wb_{ij} \le wc_{ij}$ (since $0 \le b_{ij} \le c_{ij}$).

We can distinguish between these two alternative forms by measuring the
derivative of the
concentration away from steady-state and comparing it to the derivative
predicted by each form of
the equation. This requires only a few extra thoughtfully-selected
measurements. The details are
in section S2
of the supplement.

\section*{Simulated six-gene subnetwork in mouse ESC}

To demonstrate the inference approach, we apply our method to a
synthetic six-gene system based on the Oct4, Sox2, Nanog, Cdx2, Gcnf,
Gata6 subnetwork in a mouse embryonic stem cell (ESC).
\citet{chickarmane08} developed this system based on a synthesis
of knowledge about ESC gene regulation accumulated over the past two
decades [\citet {chickarmane08}]. The network structure is shown
in Figure~\ref{fignetworks}(a), and the detailed model is given by the
following system of ODEs in the six genes:
\begin{eqnarray}
\label{eqESC}
\frac{d[O]}{dt} &=& \bigl(0.001 + [A] + 0.005 [O][S] + 0.025
[O][S][N]\bigr)\nonumber\\
&&{}/\bigl(1 + [A] + 0.001[O] + 0.005[O][S]\nonumber\\
&&\hspace*{8.2pt}{} + 0.025[O][S][N] + 10
[O][C] + 10[Gc]\bigr)\nonumber\\
&&{} - 0.1 [O],
\nonumber\\
\frac{d[S]}{dt} &=& \frac{0.001 + 0.005 [O][S] + 0.025 [O][S][N]
}{1 + 0.001[O] + 0.005[O][S] + 0.025[O][S][N]}\nonumber\\
&&{} - 0.1 [S],
\\
\frac{d[N]}{dt} &=& \frac{0.001 + 0.1[O][S] + 0.1 [O][S][N] }{1 +
0.001[O] + 0.1[O][S] + 0.1[O][S][N] + 10[O][G]}\nonumber\\
&&{} - 0.1 [N],
\nonumber\\
\frac{d[C]}{dt} &=& \frac{ 0.001 + 2[C] }{1 + 2[C] + 5[O][C]} - 0.1 [C],
\nonumber\\
\frac{d[Gc]}{dt} &=& \frac{ 0.001 + 0.1[C] + 0.1[G] }{ 1 + 0.1[C]
+ 0.1[G]} - 0.1 [Gc],
\nonumber\\
\frac{d[G]}{dt} &=& \frac{ 0.1 + [O] + 0.00025[G]}{1 + [O] + 0.00025[G]
+ 15[N]} - 0.1 [G].\nonumber
\end{eqnarray}
This model has many of the same qualitative characteristics as the
biological mouse ESC network [\citet{chickarmane08}]. In particular, the
system can support four different steady-states: embryonic stem cell
(ESC), differentiated stem cell (DSC), endoderm and trophectoderm, and
can switch from one to another when certain genes' expression levels
are changed. In the Oct4 equation, $A$~represents an external
activating factor whose concentration $[A]$ depends on the culture
condition. Each of the four steady-states has a corresponding value of
$[A]$: $10$ for ESC and DSC, $25$ for endoderm, and $1$ for
trophectoderm. For the remainder of this paper, we will regard $[A]$ as
known. The explicit system of ODEs (\ref{eqESC}) allows us to generate
data to fit our model and also to quantitatively compare our recovered
solution to the ground truth. The qualitative similarity of this
synthetic network to a real biological network gives us confidence that
our results in this numerical experiment are likely to translate well
to real biological networks.

We observe that the Cdx2, Gcnf and Gata6 equations have alternative
forms (provided we ignore the very small constant term in the $\frac
{d[C]}{dt}$ equation and $[G]$ term in the $\frac{d[G]}{dt}$). With the
minimum possible value of $w$, the alternative
forms are as follows:
\begin{eqnarray}\label{eqaltESC}
\frac{d[C]}{dt} &=& \frac{ 0.95 }{1 + 2.5[O] }\qquad (w = 2),
\nonumber\\
\frac{d[Gc]}{dt} &=& \frac{ 0.1001[Gc] + 0.01[C][Gc] +
0.01[Gc][G] }{ 1 + 0.1[Gc] + 0.01[C][Gc] + 0.01[Gc][G]}\nonumber\\
&&{} - 0.1 [Gc] \qquad(w =
0.1),
\\
\frac{d[G]}{dt} &=& \frac{0.111[G] + 0.111[O][G] }{1 + 0.111[G] +
0.111[O][G] + 1.67[N][G]} \nonumber\\
&&{}- 0.1 [G] \qquad(w = 0.111).\nonumber
\end{eqnarray}
To resolve the specific form, we will apply our method twice, once
allowing self-regulation and again disallowing it. Then we will compare
the two recovered forms of each equation and the quality of the fits to
determine whether nonidentifiability exists in each case. If so, we
will break the tie by examining derivatives.

To fit the model, we collect data on the expression levels of all six
genes at many different steady-states. First we measure the expression
levels at all four wildtype steady-states: SC, DSC, endoderm and
trophectoderm. We also induce additional perturbed steady-states by
simulating knockdowns and overexpression of each gene, based on
physical gene perturbation experiments
[\citet{rodriguez07,zafarana09}]. For a knockdown, we hold a gene
at one-fifth of its steady-state expression level; for overexpression
we hold a gene at twice its steady-state level. In each case we wait
for the system to settle to a new steady-state, then measure the
expression levels. Figure~\ref{figOKD} shows the expression
trajectories during Oct4 knockdown from the ESC steady-state as an
example.

\begin{figure}

\includegraphics{645f03.eps}

\caption{Gene expression trajectories during an Oct4 knockdown from SC
steady-state. The expression of Oct4 is artificially reduced to $20\%$
of its SC steady-state expression level and held there, causing the
expression levels of the targets of Oct4 to change in response, which
in turn impact their targets. The system eventually reaches a new
steady-state different from SC. We measure the vector of expression
levels at the new steady-state and use it as data in the inference
algorithm. Since Oct4 is knocked down, this induced steady-state does
not provide useful information about the Oct4 equation, but it is
useful for understanding the role of Oct4 and other genes in the
equations of the remaining five genes.}
\label{figOKD}
\end{figure}

\begin{figure}

\includegraphics{645f04.eps}

\caption{Recovery of a synthetic gene regulatory network based on the
biological ESC network using our inference algorithm. The diagrams
represent sytems of ODEs that quantitatively model the gene
interactions. Edge color indicates activation (green) or repression
(red), and edge weights correspond to coefficient magnitudes. The
arrows point from regulator to target, and self-loops indicate
self-regulation. The yellow star represents the third-order complex
OSN. (In addition to all possible first- and second-order terms, we
allow this special third-order term with a free coefficient.) The left
figure represents the original system of ODEs used to generate the
data. The center figure shows the network recovered using our inference
algorithm on noiseless data, and the right figure shows the recovery
with $1\%$ noise added. Both recovered networks reflect coefficient
thresholding at $0.1 \%$ (noiseless case) or $1 \%$ (noisy case) of the
largest recovered coefficent in each gene equation (with the exception
of the noiseless-case Oct4 equation, thresholded at $0.01\%$ to show
the successful recovery of weak edges). The algorithm performs almost
perfectly in the noiseless case, except for a false positive repressor
on Gata6 and two very weak activation edges missing. In the noisy case,
the algorithm recovers all of the strong edges, but misses some of the
weaker ones and returns a few small false positives at our chosen
thresholding level. Overall, the method captures the major network
structure even in the noisy case.}
\label{fignetworks}
\end{figure}

The details of the simulation are given in section S3
of the supplement. We begin by testing the algorithm on noiseless data.
We solve the optimization
problem (\ref{eqopt}) once, then we solve it again with additional
constraints prohibiting
self-regulation. In each case we use cross-validation to select the
sparsity parameter
$\lambda$ (Figure~S1).
The quality of the fit is comparable for the latter three equations
whether we restrict
self-regulation or not, while for the first three equations restricting
self-regulation has a
significant negative impact on the fit (Table S1),
indicating that the first three equations are unambiguous while the
last three have two possible
forms. To resolve the nonidentifiability in the latter three equations,
we measure the
derivatives of Cdx3, Gcnf and Gata6 immediately after some additional
informative perturbations:
Oct4, Cdx2 and Nanog knockouts, respectively (Figure S2).
The test reveals that Gcnf and Gata6 have the simple form, while Cdx2
has a higher-order form. In this example,\ the original coefficients are
recovered almost exactly:
\begin{eqnarray}
\label{eqnoiselessrecovery}
\frac{d[O]}{dt} &=& \bigl(0.001 + [A] + \bigl(0.005 [O][S] + 0.025
[O][S][N]\bigr)\bigr)\nonumber\\
&&{}/\bigl(1 + [A] + \bigl(0.001[O] + 0.005[O][S] + 0.025[O][S][N]\bigr)\nonumber\\
&&\hspace*{160.5pt}{} + 10
[O][C] + 10[Gc]\bigr)\nonumber\\
&&{} - 0.1 [O],
\nonumber\\
\frac{d[S]}{dt} &=& \frac{0.001 + 0.005 [O][S] + 0.025 [O][S][N]
}{1 + 0.005[O][S] + 0.025[O][S][N]} - 0.1 [S],
\\
\frac{d[N]}{dt} &=& \frac{0.1 [O][S] + 0.1 [O][S][N] }{1 +
0.1[O][S] + 0.1[O][S][N] + 10[O][G]} - 0.1 [N],
\nonumber\\
\frac{d[C]}{dt} &=& \frac{2[C] }{1 + 2[C] + 5[O][C]} - 0.1 [C],
\nonumber\\
\frac{d[Gc]}{dt} &=& \frac{ 0.001 + 0.1[C] + 0.1[G] }{1 + 0.1[C]
+ 0.1[G]} - 0.1 [Gc],
\nonumber\\
\frac{d[G]}{dt} &=& \frac{ 0.1 + [O]}{1 + [O] + 0.03[N][Gc] + 15[N]} - 0.1
[G].\nonumber
\end{eqnarray}
Next we add zero-mean Gaussian noise to each measurement, with standard
deviation $1 \%$ of the
measurement magnitude. We use the same steady-states as in the
noiseless case, plus
overexpression-knockdown of each pair of genes starting from ESC and
DSC. Using a similar approach
(detailed in section S3
of the supplement), we recover:
\begin{eqnarray} \label{eqnoisyrecovery}
\frac{d[O]}{dt} &=& \frac{ [A] }{1 + [A] + 9.9[Gc] + 9.9 [O][C] } - 0.1
[O],
\nonumber\\[-2pt]
\frac{d[S]}{dt} &=& \frac{0.001 [O][S] + 0.0005[S][N] + 0.025
[O][S][N] }{1 + 0.001[O][S] + 0.0005[S][N] + 0.025[O][S][N]} - 0.1 [S],
\nonumber\\[-2pt]
\frac{d[N]}{dt} &=& \frac{0.09[O][S][N] }{1 + 0.1[G][Gc] +
0.09[O][S][N] + 9.1[O][G]} - 0.1 [N],
\nonumber\\[-9pt]\\[-9pt]
\frac{d[C]}{dt} &=& \frac{ 2[C] }{1 + 2[C] + 5[O][C] } - 0.1 [C],
\nonumber\\[-2pt]
\qquad\frac{d[Gc]}{dt} &=& \frac{ 0.1[C] + 0.1[G] }{ 1 + 0.1[C] +
0.1[G]} - 0.1 [Gc],
\nonumber\\[-2pt]
\frac{d[G]}{dt} &=& \frac{ 0.1 + 0.9[O] }{1 + 0.9[O] + 14.2[N] } - 0.1
[G].\nonumber
\end{eqnarray}
In order to produce clean equations and network diagrams, we choose
appropriate thresholds for each equation below which we zero the
coefficients. (In practice, choosing thresholds is a judgment call
based on the expected number of regulators, the noise level of the data
and the level of detail appropriate for the application.) We set the
thresholds at $0.1\%$ (noiseless case) or $1\%$ (noisy case) of the
largest coefficient recovered for each equation. For example, the
largest recovered coefficient in the $d[G]/dt$ equation is roughly 15
in either case, so we zero the coefficients that fall below $0.015$
(noiseless case) or $0.15$ (noisy case). The recovered systems of
equations shown above reflect these choices. In the noiseless case,
relaxing the threshold on the Oct4 equation to $0.01\%$ leads to the
recovery of more correct terms, listed in parentheses. For
completeness, we also provide receiver operating characteristic (ROC)
curves in Figure~\ref{figROC} to show the trade-off between true
positives and false positives at other thresholds. The network diagrams
in Figure~\ref{fignetworks}(b), (c) include an edge if the corresponding
\begin{figure}

\includegraphics{645f05.eps}

\caption{ROC curves for recovered networks from noiseless (left) and
noisy (right) data showing the trade-off between true positive rate
(TPR) and false positive rate (FPR) for edge recovery. The ROC curves
show the TPR and FPR that result from a range of coefficient threshold
choices above which we consider an edge to have been recovered. For the
equation $dx_i/dt = \cdots$ and threshold $t$, TPR is defined as the
proportion of true edges $j$ with $c_{ij}^{\mathrm{recovered}} > t$ and FPR
as the proportion of false edges with $c_{ij}^{\mathrm{recovered}} >t$. For
equations with two possible forms, we compare the simple forms of the
true and recovered equations. Each gene equation has a different ROC
curve as indicated by the legend. The dotted black line is the expected
ROC curve for the ``random guessing'' algorithm, while the $(0,1)$ point
corresponds to a perfect algorithm (in fact, our algorithm performs
perfectly for the Gcnf equation).}
\label{figROC}
\end{figure}
coefficient is above the threshold, with weights reflecting the size of
the coefficients. These diagrams show that the recovery is nearly
perfect in the noiseless case: using the gentler threshold for the Oct4
equation, we recover all the true edges except for three very weak
ones, and return just one small false positive repressor in the Gata6
equation. In the noisy case, we recover all the large coefficients
correctly, although there are a few small false positives and we miss
several of the weakest edges. Overall, the method is able to capture
the major network structure.

\section*{Discussion}

Our experiment on the synthetic ESC system demonstrates that our
algorithm can be used to infer a complex dynamical systems model of
gene regulation and that the method can tolerate low levels of noise.
Term selection from among all possible single gene and gene-complex
regulators (up to second-degree interactions, plus the third-degree
interaction OSN) was successful. The inferred equations are easy to
interpret in terms of gene networks, and the detailed quantitative
information allows for prediction of future expression trajectories
from any starting point.\vadjust{\goodbreak}

The approach is also scalable. Since we have formulated our problem as
a convex optimization problem [equation (\ref{eqopt})], it can be solved
efficiently even for large systems using prepackaged software.
Furthermore, it is trivially parallelizable, since we need to solve a
version of (\ref{eqopt}) to infer the differential equation
coefficients $b_{ij}, c_{ij}$ for each gene $i$. Parallelization is
even more helpful for the cross-validation step, where we need to solve
equation (\ref{eqopt}) for each gene and a sequence of choices sparsity
parameter $\lambda$. We tested the scalability by running the algorithm
with the parallelization discussed above on a simulated 100-gene
system. The algorithm ran correctly in a reasonable time frame (a few
hours) on a computing cluster.

The high resolution of our model is one of its most valuable features,
but it means that accurate term selection may require much data,
especially in the presence of noise. In our experiment, when we added
$1\%$ Gaussian noise, we needed extra data (knockdown/overexpression
pairs) in order to accurately select terms. When we tried $5\%$ noise,
the algorithm consistently selected the large terms in five of the six
equations, but we had to add even more data in order to correctly
identify the major repressor in the Nanog equation. The Nanog equation
is subtle in that Oct4 acts as both an activator in complexes with Sox2
and Nanog and a repressor in a complex with Gata6, so the algorithm
tends to select different Gata6 complexes (or the Gata6 singleton) as
the major repressor when the data is insufficient. In the $5\%$ noise
case, we needed additional data on the role of Gata6 (double-knockdowns
and double-overexpression of pairs including Gata6 from ESC and DSC) in
order to select Oct4-Gata6 as the major repressor of Nanog fairly
consistently. As discussed earlier, another difficulty is the
nonidentifiability that arises from accounting for self-regulation
while restricting data to steady-states. Distinguishing between the two
possible forms of nonidentifiable equations requires extra derivative
data (which can be collected experimentally, although it is more
difficult and time-consuming) and extra steps in the algorithm. The
constraints on the convex optimization problem (\ref{eqopt}), which
arise from thermodynamic considerations, are sufficient to prevent
further nonidentifiability, but in certain cases, certain problems can
suffer from near-nonidentifiability of other forms, which may
contribute to the challenge of term-selection with noisy or limited
data. We ensure accurate term selection by making sure we include
enough diverse, high-quality steady-state measurements.

We should also note that our model does not account for the intrinsic
noise in gene transcription and translation, although these processes
are inherently stochastic, since TF and RNAP binding result from chance
collisions between molecules in the cell. However, the stochastic
version of our rational-form transcription model is highly complex and
there is currently no satisfactory method for its inference. Studying
the deterministic evolution plus additive noise is standard practice
for all but linear models of gene expression, and treatment of the
deterministic model provides insight into the stochastic model. Here we
focus on the additive noise case and leave the study of intrinsic noise
for future investigation.

\section*{Conclusions}

The model we use is based on the detailed thermodynamics of gene
transcription, and quantitatively captures the full spectrum of
regulatory phenomena in a detailed, physically interpretable,
predictive manner. Since we can formulate the model fitting problem as
a convex optimization problem, we can solve it efficiently and scalably
using prepackaged software. $\ell_1$-regularization allows for
term-selection while maintaining the problem convexity. The experiments
required to collect the necessary steady-state gene expression data are
straightforward to perform, as technologies for knockdowns and
overexpression are well-established and measuring gene expression is
relatively simple. The model accounts for activation and repression by
single-protein TFs and synergistic complexes as well as
self-regulation, and describes the magnitude of each type of regulation
in quantitative detail. Furthermore, the model can be extended to
account for environmental effects and auxiliary proteins involved in
regulation, including enhancers and chromatin remodelers. The fitted
model can predict the evolution of the system from any starting point.
Given a set of steady-states gene expression measurements, our
algorithm can be used to fit a model which not only predicts further
steady-states of the system, but also fully describes the transitions
between them. Finally, beside the study of gene regulation, our
approach will be useful in many other application areas where it is
necessary to infer a nonlinear dynamical system by suitable
experimentation and statistical analysis.

\begin{appendix}\label{app}
\section*{Appendix: Thermodynamic model}

In (\ref{eqexpression}), the function $f_i(y)$ represents the
probability that RNAP binds to the $i$th gene promoter. We claim that
$f_i(y)$ has the form
\[
f_i(y) \equiv p_{\mathrm{bound}}^{(i)}(y) =
\frac{ \sum_j e^{-\beta\Delta
\varepsilon_{ij}^{\mathrm{RNAP}}} P e^{-\beta\Delta\varepsilon_{ij}} \Pi_{k \in
S_{ij}} y_k} {
\sum_j (1+e^{-\beta\Delta\varepsilon_{ij}^{\mathrm{RNAP}}} P ) e^{-\beta
\Delta\varepsilon_{ij}} \Pi_{k \in S_{ij}} y_k},
\]
where $\Delta\varepsilon_{ij}$ is the binding energy of the $j$th complex
to the promoter, $\Delta\varepsilon_{ij}^{\mathrm{RNAP}}$ is the binding
energy of RNAP to the $j$th promoter-bound complex, and $P, x_j$ are
the concentrations of RNAP and gene product $j$ [Bintu et~al.
(\citeyear{bintu05apps,bintu05models})].

Any type of regulator (including no regulator at all) can be
represented in this framework. For no regulator, we take $S_{ij} =
\varnothing$ with the convention that $\Pi_{k \in\varnothing} y_k = 1$,
set $\Delta\varepsilon_{ij} =0$, and take $\Delta\varepsilon_{ij}^
{\mathrm{RNAP}}$ as the base binding energy of RNAP to the promoter. For a
repressor, $\Delta\varepsilon_{ij} < 0$ and $\Delta\varepsilon_{ij}^
{\mathrm{RNAP}} > 0$; for an activator, $\Delta\varepsilon_{ij} < 0$ and $\Delta
\varepsilon_{ij}^{\mathrm{RNAP}} < 0$.

Setting
\begin{eqnarray*}
b_{ij} &=& e^{-\beta\Delta\varepsilon_{ij}^{\mathrm{RNAP}}} P e^{-\beta\Delta
\varepsilon_{ij}},
\\
c_{ij} &=& \bigl(1+e^{-\beta\Delta\varepsilon_{ij}^{\mathrm{RNAP}}} P \bigr) e^{-\beta
\Delta\varepsilon_{ij}},
\end{eqnarray*}
we obtain the form given in Section~1:
\[
f_i(y) = \frac{ b_{ij} \Pi_{k \in S_{ij}} y_k } { \sum_j c_{ij} \Pi_{k
\in S_{ij}} y_k }.
\]
Constant terms in the numerator and denominator correspond to the
no-regulator case. Letting $c_{i0}$ denote the constant appearing in
the denominator, our convention will be to divide all of the
coefficients in the numerator and denominator by $c_{i0}$ so that the
constant 1 appears in the denominator.

\subsection*{Simplified derivation}
The derivation we present here follows Bintu et al. and Garcia et al.
[Bintu et~al. (\citeyear{bintu05apps,bintu05models}),
\citet{garcia11}]. For simplicity, we will prove the following
claim for the simplified case with one regulator $y_1$ (as well as the
possibility of RNAP binding with no regulator):
\[
p_{\mathrm{bound}}^{(i)} = \frac{ e^{-\beta\Delta\varepsilon_{i0}^
{\mathrm{RNAP}}} p + e^{-\beta\Delta\varepsilon_{i1}^{\mathrm{RNAP}}} p e^{-\beta
\Delta\varepsilon_{i1}} y_1 } {
(1 + e^{-\beta\Delta\varepsilon_{i0}^{\mathrm{RNAP}}} p) + (1+e^{-\beta
\Delta\varepsilon_{ij}^{\mathrm{RNAP}}} p ) e^{-\beta\Delta\varepsilon_{i1}}
y_1 }.\vadjust{\goodbreak}
\]

We will use the following notation:
$\varepsilon_{P,i1}^S$ is the energy of the state in which RNAP is
specifically bound to the regulator-promoter complex, $\varepsilon
_{P,i0}^S$ is the energy of the state in which RNAP is specifically
bound to the promoter without the regulator,
$\varepsilon_{P}^{\mathrm{NS}}$ is the energy when RNAP is bound to a nonspecific
binding site, $\varepsilon_{i1}^S$ is the energy when $y_1$ is
specifically bound to the promoter, and $\varepsilon_{i1}^{\mathrm{NS}}$ is energy
when $y_1$ is bound to a nonspecific binding site. Then
\begin{eqnarray*}
\Delta\varepsilon_{i0}^{\mathrm{RNAP}} &=& \Delta\varepsilon_{P,i0}
\equiv \varepsilon_{P,i0}^S - \varepsilon_{P}^{\mathrm{NS}},\\
\Delta\varepsilon_{i1}^{\mathrm{RNAP}} &=& \Delta\varepsilon_{P,i1}
\equiv \varepsilon_{P,i1}^S - \varepsilon_{P}^{\mathrm{NS}},\qquad
\Delta\varepsilon_{i1} \equiv\varepsilon_{y_1}^{S} -
\varepsilon_{y_1}^{\mathrm{NS}}.
\end{eqnarray*}
Suppose that we have $j$ RNA polymerase molecules and $k$ molecules of
gene product $1$ (the regulator). We model the genome as a
``reservoir''
with $n$ nonspecific binding sites (to which either RNAP or regulator
can bind). One of these sites is the promoter of gene $i$. Four
different classes of configurations interest us:
\begin{enumerate}
\item empty promoter,
\item regulator bound to promoter,
\item regulator and RNAP bound to promoter,
\item RNAP only bound to promoter.
\end{enumerate}
These correspond to the following partial partition functions, which
represent the ``unnormalized probabilities'' of each configuration:
\begin{enumerate}
\item$Z(j,k)$,
\item$Z(j,k-1)e^{-\beta\varepsilon_{i1}^S}$,
\item$Z(j-1,k-1)e^{-\beta\varepsilon_{i1}^S}
e^{-\beta\varepsilon_{P,i1}^S}$,
\item$Z(j-1,k) e^{-\beta\varepsilon_{P,i0}^S}$,
\end{enumerate}
where $Z(j,k) = \frac{n!}{j!k!(n-j-k)!} e^{-\beta r \varepsilon_{i1}^{\mathrm{NS}}}
e^{-\beta\varepsilon_{P}^{\mathrm{NS}}}$.

$Z(j,k)$ is equal to the total number of arragements of RNAP and
regulator on the nonspecific binding sites times the Boltzmann factor,
which gives the relative probability $e^{-\beta\varepsilon}$ of a
particular state in terms of its energy $\varepsilon$.

Since RNAP binds the promoter only in the third and fourth classes of
configurations, the probability that RNAP binds the promoter is equal
to the unnormalized probability of the third and fourth configurations
divided by the ``total probability'' (the sum of the unnormalized
probabilities of all classes of configurations). Hence,
\begin{eqnarray*}
p_{\mathrm{bound}} & = & \bigl(Z(j-1,k) e^{-\beta\varepsilon_{P,i0}^S} + Z(j-1,k-1)e^{-\beta
\varepsilon_{i1}^S} e^{-\beta\varepsilon_{P,i1}^S}\bigr)\\
&&{}/ \bigl(
Z(j,k) + Z(j-1,k) e^{-\beta\varepsilon_{P,i0}^S}\\
&&\hspace*{6.8pt}{}  + Z(j,k-1)e^{-\beta
\varepsilon_{i1}^S}+ Z(j-1,k-1)e^{-\beta\varepsilon_{i1}^S} e^{-\beta
\varepsilon_{P,i1}^S} \bigr)
\\
&\approx&\biggl(\frac{n^{j-1}n^{k}}{(j-1)!k!} e^{-\beta k \varepsilon
_{i1}^{\mathrm{NS}}} e^{-\beta(j-1) \varepsilon_{P}^{\mathrm{NS}}}e^{-\beta\varepsilon
_{P,i0}^S}\\
&&\hspace*{5pt}{} + \frac{n^{j-1}n^{k-1}}{(j-1)!(k-1)!} e^{-\beta(k-1)
\varepsilon_{i1}^{\mathrm{NS}}} e^{-\beta(j-1) \varepsilon_{P}^{\mathrm{NS}}} e^{-\beta
\varepsilon_{i1}^S} e^{-\beta\varepsilon_{P,i1}^S}\biggr)\\
&&{}\bigg/
\biggl(\frac{n^{j}n^{k}}{j!k!} e^{-\beta k \varepsilon_{i1}^{\mathrm{NS}}} e^{-\beta j
\varepsilon_{P}^{\mathrm{NS}}}\\
&&\hspace*{16pt}{}
+ \frac{n^{j-1}n^{k}}{(j-1)!k!} e^{-\beta k \varepsilon_{i1}^{\mathrm{NS}}}
e^{-\beta(j-1) \varepsilon_{P}^{\mathrm{NS}}}e^{-\beta\varepsilon_{P,i0}^S} + \cdots\biggr)
\\
& =&
\biggl(\frac{j}{n}e^{\beta\varepsilon_{P}^{\mathrm{NS}}}e^{-\beta\varepsilon
_{P,i0}^S} + \frac{j}{n} \frac{k}{n} e^{\beta\varepsilon_{i1}^{\mathrm{NS}}}
e^{\beta\varepsilon_{P}^{\mathrm{NS}}} e^{-\beta\varepsilon_{i1}^S} e^{-\beta
\varepsilon_{P,i1}^S}\biggr)\\
&&{}\bigg/\biggl(1 + \frac{j}{n}e^{\beta\varepsilon_{P}^{\mathrm{NS}}}e^{-\beta\varepsilon
_{P,i0}^S} + \frac{k}{n} e^{\beta\varepsilon_{i1}^{\mathrm{NS}}}
e^{-\beta\varepsilon_{i1}^S}\\
&&\hspace*{37pt}{}+ \frac{j}{n} \frac{k}{n} e^{\beta\varepsilon_{i1}^{\mathrm{NS}}} e^{\beta\varepsilon
_{P}^{\mathrm{NS}}} e^{-\beta\varepsilon_{i1}^S} e^{-\beta\varepsilon_{P,i1}^S}\biggr)
\\
& = & \frac{ ({j}/{n}) e^{-\beta\Delta\varepsilon_{P,i0}} + ({j}/{n})
({k}/{n}) e^{ -\beta\Delta\varepsilon_{i1}} e^{-\beta\Delta\varepsilon
_{P,i1}} } {
1 + ({j}/{n}) e^{-\beta\Delta\varepsilon_{P,i0}}
+ ({k}/{n}) e^{ -\beta\Delta\varepsilon_{i1}}
+ ({j}/{n}) ({k}/{n}) e^{ -\beta\Delta\varepsilon_{i1}} e^{-\beta
\Delta\varepsilon_{P,i1}} }
\\
& = & \frac{ ({j}/{n}) e^{-\beta\Delta\varepsilon_{P,i0}} + ({j}/{n})
({k}/{n}) e^{ -\beta\Delta\varepsilon_{i1}} e^{-\beta\Delta\varepsilon
_{P,i1}} } {
1 + ({j}/{n}) e^{-\beta\Delta\varepsilon_{P,i0}} + ({k}/{n}) e^{
-\beta\Delta\varepsilon_{i1}} ( 1 + ({j}/{n}) e^{-\beta\Delta\varepsilon
_{P,i1}} ) }
\\
& = & \frac{ p e^{-\beta\Delta\varepsilon_{i0}^{\mathrm{RNAP}}} + p y_1 e^{
-\beta\Delta\varepsilon_{i1}} e^{-\beta\Delta\varepsilon_{i1}^
{\mathrm{RNAP}}} } {
1 + p e^{-\beta\Delta\varepsilon_{i0}^{\mathrm{RNAP}}} + y_1 e^{ -\beta
\Delta\varepsilon_{i1}} ( 1 + p e^{-\beta\Delta\varepsilon_{i1}^
{\mathrm{RNAP}}} ) },
\end{eqnarray*}
where in the second line we used the approximation $\frac
{n!}{j!k!(n-j-k)!} \approx\frac{n^{j}n^{k}}{j!k!} $ which holds for
$j, k \ll n$, in the third we divided by $\frac{n^{j}n^{k}}{j!k!}
e^{-\beta k \varepsilon_{i1}^{\mathrm{NS}}}e^{-\beta\varepsilon_{P}^{\mathrm{NS}}}$, in the
fourth we used the identities $\Delta\varepsilon_{P,i0} = \varepsilon
_{P,i0}^S - \varepsilon_{P}^{\mathrm{NS}}$, $\Delta\varepsilon_{P,i1} = \varepsilon
_{P,i1}^S - \varepsilon_{P}^{\mathrm{NS}}$, $\Delta\varepsilon_{i1} \equiv\varepsilon
_{i1}^{S} - \varepsilon_{i1}^{\mathrm{NS}}$, and in the last we substituted in the
definitions $\frac{j}{n} = p$, $\frac{k}{n} = y_1$, $\Delta\varepsilon
_{i0}^{\mathrm{RNAP}} = \Delta\varepsilon_{P,i0}$, $\Delta\varepsilon_{i1}^
{\mathrm{RNAP}} = \Delta\varepsilon_{P,i1}$.
\end{appendix}

\section*{Acknowledgments}

Many thanks to Xi Chen for helpful discussions.

\begin{supplement
\stitle{Nonidentifiability, tie-breaking and synthetic network study details}
\slink[doi]{10.1214/13-AOAS645SUPP}
\sdatatype{.pdf}
\sfilename{aoas465\_supp.pdf}
\sdescription{We discuss nonidentifiability and tie-breaking in
Sections S1
and S2
by proving the equivalence of two different equation forms at
steady-state and describing methods
for determining the true form of an ambiguous equation. In Section S3
we provide the details of our study of a simulated six-gene network in
mouse ESC, including parameter selection, tie-breaking and
thresholding.}
\end{supplement}


\section{Introduction}
\label{sec:intro}
The origin of ultra-high-energy (UHE) cosmic rays (CRs), that is of cosmic
particles with energies $E \gtrsim 10^{19}$~eV, is presently unknown.
However, there are numerous hints, both in data and in theory, which might
help to constrain possible models of UHECR sources. In particular, the
observation \cite{HiRes-GZK, PAO-GZK, TA-GZK} of the spectrum steepening
consistent with the Greisen--Zatsepin--Kuzmin \cite{G, ZK} (GZK) cutoff,
together with the global isotropy of the arrival directions
(see e.g.\ ~\cite{TA-anisotropy}) and the fact that the UHE
particles are not expected to be confined by the Milky-Way magnetic field
(see e.g.\ \cite{NotConfined}) suggest that the bulk of cosmic rays at
these energies have extragalactic origin\footnote{See however Ref.\
\cite{STdoublet} and references therein where possible exceptions are
discussed.}. Next, the lack of clustering of arrival directions at small
scales is a powerful tool \cite{DubovskyClustering} to constrain the
number density of sources. In recent data of the Pierre Auger Observatory
(PAO), no evidence for clustering at $E \gtrsim 5 \times 10^{19}$~eV is
seen \cite{PAO-clustering} which translates into the number density of $n
\gtrsim 10^{-4}$ sources per cubic Megaparsec, which means that sources of
these extreme particles should not be exceptional and, most probably, some
of them should be located relatively nearby. The latter fact gets further,
though limited by statistical significance, support from the shape of the
GZK feature in the spectrum which does not seem to be very sharp (cf.\
Ref.~\cite{SemikozMild}).

On the other hand, it is a nontrivial task to find particular
astrophysical objects which could serve as UHECR accelerators. Even
without a detailed modelling of the acceleration process, a number of
simple estimates rule out many classes of potential sources. These
simple criteria include in particular the geometrical (Hillas) criterion
\cite{Hillas} and estimates of radiative energy losses of particles being
accelerated (see e.g.\ Refs.~\cite{Aha-losses, Medv-losses}). Analysis of
the modern astrophysical data demonstrates \cite{UFN-HillasPlot} that the
combination of these constraints leaves just a few candidate classes of
sources capable of acceleration of particles to UHE energies. Leaving
aside large-scale structures where interaction losses are expected to
suppress the energy gain, the conventional diffusive (e.g., relativistic
or non-relativistic shock) acceleration may work only in ultrarelativistic
jets, hot spots and lobes of exceptional active galaxies (powerful radio
galaxies and blazars) which are not that abundant in the nearby Universe.
For very special field configurations when synchrotron losses are
suppressed and the curvature radiation dominates, possible acceleration
sites include also gamma-ray bursts (GRBs) and immediate neighbourhood of
supermassive black holes (SMBHs) in the galactic nuclei. While it is
unclear whether these field configurations may be present in GRBs, recent
IceCube results disfavour the GRB scenario anyway \cite{IceCubeGRB} (see
however Ref.~\cite{Dar}). At this level of reasoning, the SMBH environment
remains a viable option.

A natural assumption is that numerous UHECR sources are not
identical -- there should be less and more powerful accelerators where the
maximal energies, injection spectra and fluxes of accelerated particles are
different. Until now, numerous attempts to model the sources of UHECRs and
to confront theoretical predictions with experiments often assumed that
these parameters are fixed once and for the entire Universe (see e.g.\
Refs.~\cite{Berezinsky, Oleg} and numerous other works; see however
Ref.~\cite{Ptuskin} where acceleration in non-identical jets was
considered). While, for numerous sources, the assumption of equal fluxes is
well justifiable (in the sence that only the mean flux of a large sample
of sources is important and this mean flux does not vary significantly
from one region in the Universe to another) and the injection spectrum is
often fixed by the acceleration model, the maximal energies are expected
to vary significantly. As it was recognised in Ref.~\cite{SemikozMaxEn},
these variations affect the observable spectrum seriously. In this work,
we attempt to present a toy model of numerous and different sources of
UHECRs and, within certain assumptions, to confront it with the
experimental data.

To this end, we choose a simple toy model of particle acceleration in the
immediate vicinity of SMBH put forward in Refs.~\cite{NerSem, NerSemTk}.
The reason to choose this particular model is twofold. First, unlike many
other models, it allows~\cite{UFN-HillasPlot} for UHECR acceleration in
numerous nearby sources. Second, as we will see below, within some
realistic assumptions, the acceleration capabilities of a source are
determined by a single parameter, the SMBH mass. At the same time, the
demography of SMBHs is well studied by astrophysicists and we take this
advantage to describe the population of sources easily.

The rest of the paper is organized as follows. In Sec.~\ref{sec:BH}, we
give a brief review of the acceleration model of Refs.~\cite{NerSem,
NerSemTk} and make a bridge between the parameters which determine the
maximal energy of accelerated particles and the SMBH mass. In
Sec.\ref{sec:population}, we discuss the astrophysical data on the SMBH
population, merge them with the acceleration model, calculate the spectrum
of UHECRs with the account of propagation from source to the observer and
compare it with the experimental cosmic-ray data. We obtain a good
agreement with the observed spectrum by fitting the spectrum with only
two continuous parameters, the overall normalization and a single free
parameter of the model. Sec.~\ref{sec:constraints} demonstrates that the
population model of Sec.~\ref{sec:population} satisfies simple
observational constraints: it does not produce too much secondary gamma
rays, it results in an acceptable number density of sources and, with the
best-fit normalization, it does not require enormous luminosity of a
single source. We give our conclusions and discuss our results in
Sec.~\ref{sec:concl}.

\section{A toy model of particle acceleration in the black-hole
magnetosphere}
\label{sec:BH}
A toy model of particle acceleration in the black-hole magnetosphere was
proposed by Neronov et al.~\cite{NerSem, NerSemTk}. Let us briefly discuss
the model and its parameters.

Assume that a stationary rotating black hole without
electric charge is embedded into the external magnetic field, homogeneous
at the horizon distance scale. In general, the magnetic field is inclined
at some angle $\chi$ with respect to the black-hole rotation axis. There is
a well-known exact solution of Maxwell's equations in the Kerr metric for
each inclination angle $\chi$ of an asymptotically homogeneous magnetic
field \cite{BHsolution1, BHsolution2}. For instance, if $\chi=0$, then a
rotation-induced electric field is parallel to the magnetic one on the
symmetry axis and its direction depends on the directions of both the
magnetic field and the black hole's rotation velocity. Thus in the region
near the rotation axis, particles moving along magnetic lines are
accelerated by the electric field.

In this case, radial components of the electric and magnetic fields in
units $\hbar=c=G=1$ in locally non-rotating frame in Boyer-Lindquist
coordinates on the symmetry axis are:
$$
B_{\widehat{r}}=B_{0}\left(1-\frac{4a^2Mr}{(r^2+a^2)^2}\right),
$$
\begin{equation}
\label{eq:1}
E_{\widehat{r}}=-\frac{2aMB_{0}(r^2-a^2)}{(r^2+a^2)^2},
\end{equation}
where $M$ is the black hole's mass, $a \leq M$  is its angular momentum
per unit mass, $r$ is the radial coordinate and $B_{0}$ is the
external homogeneous magnetic field.

Neglecting for the moment the energy losses, the maximal energy gain of the
accelerating particle with a charge $Ze$ is determined by the available
potential difference in the region along the rotation axis,
\begin{equation}
\label{eq:2}
{\cal E}_{\rm max}(a)=\int_{r_{\rm hor}}^{r_{\rm max}} Ze E(r,a)dr,
\end{equation}
$r_{\rm hor}=M+\sqrt{M^2-a^2}$ is the radius of the black-hole horizon
where particle acceleration starts while $r_{\rm max}$ limits the
size of the region along the rotation axis, where acceleration is possible
(in Ref.~\cite{NerSemTk}, it is called ``the vacuum gap'' due to the
absence of numerous charged particles in this region except of the
single test particle being accelerated, which does not change
the electromagnetic field).

Of course, one can rewrite Eq.~(\ref{eq:2}) for the potential difference
and ${\cal E}_{\rm max}$ in terms of the distance-averaged electric field
$\bar{E}$,
$$
{\cal E}_{\rm max}(a)= Ze \bar{E}(a) H ,
$$
$$\bar{E}(a)=\frac{1}{H}\int_{r_{\rm hor}}^{r_{\rm max}} Ze E(r,a)dr,$$
where $H={r_{\rm hor}}-{r_{\rm max}}$.

By making use of  Eq.~(\ref{eq:1}) one obtains
$$
{\cal E}_{\rm max}(a)=2ZeMB_{0}\left.\frac{ra}{r^2+a^2}\right|_{r_{\rm hor}}^{r_{\rm
max}}
$$
\begin{equation}
\label{eq:4}
=2aZeMB_{0} \left[
\frac{(r_{\rm hor}+H)}{(r_{\rm hor}+H)^2+a^2}-\frac{r_{\rm hor}}{r_{\rm
hor}^2+a^2}\right].
\end{equation}
It is easy to see that the expression in square brackets equals to
$-H(H+2Mr_{\rm hor}-2a^2)$ and therefore is always negative
for any value of $H$ (remember that the angular momentum per unit mass $a$
varies between zero and $M$).
Thus for parallel magnetic field and angular momentum, $aM>0$ and the
radial component of the electric field $E_{\widehat{r}}$ on the black-hole
rotation axis is negative. So if there are negative charges near the
rotation axis, they will be accelerated away from the black hole.
According to Eq.~(\ref{eq:4}) the difference between their energies at
$r_{\rm max}$ and $r_{hor}$ is positive, so they gain energy while moving
along the rotation axis away from the black hole.
In the opposite case, when the magnetic field and angular momentum are
antiparallel, the radial component of the electric field on the black hole
rotation axis is positive.  So positive-charge particles, situated near the
rotation axis, are accelerated away from the black hole; according to
Eq.~(\ref{eq:4}), the energy gain in this
case is positive for the positive-charge particles.

In Ref.~\cite{NerSemTk},
a simple expression ${\cal E}_{\rm max}\sim \  ZeB_{0}H$
rather than Eq.~(\ref{eq:4}) was used; however,
more precisely,
Eq.~(\ref{eq:4}) implies ${\cal E}_{\rm max}< \  ZeB_{0}H$.
The dependence of $\xi={\cal E}_{\rm max}/ZeB_{0}H$ from the angular
momentum $a$ for different $H$ is shown in Fig.~\ref{f_1}.
\begin{figure}
\centering \includegraphics[width=0.48\textwidth]{potential.eps}
\caption{\label{f_1}
Dependence of the normalized maximal energy $\xi$ from the angular
momentum $a$ of the black hole with mass $M$ for different values of the
vacuum-gap size $H$ (curves labeled by values of $H/M$; see the text for
definitions of parameters). }
\end{figure}
For black holes with angular momentum $a>0.1M$ and $H\sim
(1-2)R_{\rm S}=(2-4)M$, we have
${\cal E}_{\rm max}\sim 0.1\  ZeB_{0}H$; however, for slowly rotating
black holes, $a<0.1M$, $\ {\cal E}_{\rm max}$ varies from zero to
$\ \sim 0.01 ZeB_{0}H$. Hereafter, $R_{\rm S}$ denotes the Schwarzschild
radius.

However, the precise value of ${\cal E}_{\rm max}$ is often
irrelevant since particles cannot achieve this maximal energy because of
inevitable energy losses associated with the accelerated motion of the
particle. The particle energy is determined by the balance between the
energy losses and the energy gain per unit time,
$$\frac{d{\cal E}_{+}}{dt}=\frac{d{\cal E}_{-}}{dt}.$$
It was shown in Ref.~\cite{NerSemTk} that protons can be accelerated to
the energies of about $10^{20}$~eV only if the magnetic field is almost
aligned with the rotation axis.
In this case, only the curvature radiation is relevant for an accelerated
particle, and not the synchrotron one. We should also note here that we
do not consider energy losses related to interactions of accelerated
particles in the source. The maximal energy of an accelerated particle
is
$${\cal E}_{\rm curv}=\left(\frac{3}{2}\right)^{1/4}\frac{A}{Z^{1/4}}\
\frac{m}{e^{1/4}}\ E^{1/4} \ R^{1/2},$$
\begin{equation}
{\cal E}_{\rm curv} \approx 1.23 \times 10^{22}~\mathrm{eV} \ \frac{A}{Z^{1/4}}\
\left(\frac{B_0}{1~\mathrm{G}}\right)^{1/4}
\left(\frac{R}{1~\mathrm{kpc}}\right)^{1/2}{{\kappa}^{1/4}},
\label{Eq:Ecurv}
\end{equation}
where $B_{0}$ is the external magnetic field, $\ R\sim R_{\rm S}/\chi$  is
the curvature radius of magnetic-field lines, $Ze$ is the particle charge,
$Am$ is the particle mass ($A$ is the atomic number and $m$ is the nucleon
mass) and $\kappa$ is a coefficient between the electric field and
the external magnetic field, $|E_{\widehat{r}}|=\kappa B_{0}$ (see
Eq.~(\ref{eq:1})). Note that $\kappa$ is a function of the angular
momentum $a$ and the coordinate $r$. We point out that numerically,
${\cal E}_{\rm curv}<{\cal E}_{\rm max}$, cf.\ Ref.~\cite{UFN-HillasPlot},
and therefore it is ${\cal E}_{\rm curv}$ which determines the actual
maximal energy. For simplicity, we assume that all of the particles start
with equal initial conditions and so all of them are accelerated to the
same energy ${\cal E}_{\rm curv}$.

The original model of Ref.~\cite{NerSemTk} treated the magnetic field
$B_{0}$ as a free parameter. However, one may note that the field is
constrained and, in particular, cannot be too high (see
Ref.~\cite{UFN-HillasPlot} for a detailed discussion).
The maximal value of the magnetic field is determined \cite{Znajeck,
Abramovich} by the so-called Eddington limit,
$B_{\rm Ed}=10^{4}\left(\frac{M}{10^9 M_{\bigodot}}\right)^{-1/2}
\text{G}$. Though this estimate may be obtained in several different ways,
its most transparent meaning is that the magnetic-field energy density is
equal to that of the accreting plasma, corresponding to the Eddington
luminosity. To obtain the maximal possible particle energy, we should
assume that the external magnetic field is equal to the Eddington limit,
$B_{0}=B_{\rm Ed}$.
However, below, we will find the spectrum of cosmic rays in the frameworks
of this model. For this purpose we are interested in the actual maximal
particle energy as a function of black-hole mass $M$ rather than in the
upper limit. So we should recognize that the realistic magnetic field can
differ from the Eddington limit. In a general case, we can parametrize the
external magnetic field as follows:
$$B_{0}=k B_{\rm Ed}
\left(
\frac{M}{10^{9}M_{\bigodot}}
\right)^{\alpha},$$
where $\alpha$ and $k$ are some parameters. We note in passing that
several realistic models predict this kind of dependence, e.g.\ the
Shakura--Syunyaev model \cite{ShakuraSyunyaev, NovikovThorne} ($k\approx
0.31$, $\alpha=0$) or the model of Ref.~\cite{OtherBmodel} ($k\approx
0.0093$, $\alpha \approx -0.31$; see Fig.~2 of Ref.~\cite{UFN-HillasPlot}
for comparison with scarce observational data). We will consider these
options below.
We have
\begin{equation}
{\cal E}_{curv} \approx 2.9\times 10^{20}~\mathrm{eV} \ \frac{A}{Z^{1/4}}\!
\left(\frac{M}{10^{9}M_{\bigodot}}\right)^{\frac{3}{8}+\frac{\alpha}{4}}\!
{\left(\frac{\chi}{1^{\circ}}\right)}^{-\frac{1}{2}}\!(k\kappa)^{\frac{1}{4}}.
\label{Eq:EfromMass}
\end{equation}

During acceleration, the particle emits curvature photons. In what
follows, we will need to obtain an upper bound on this emission. The peak
energy of the photons is determined by the particle energy ${\cal E}$,
$$
{\cal E}_{\gamma}=\frac{3}{2}\frac{{\cal E}^3}{m^3 R},
$$
and for the upper limit, we take ${\cal E}={\cal E}_{\rm curv}$,
$$
{\cal E_{\gamma}}\sim 14~\mbox{TeV}~\frac{A^3}{Z^{3/4}}\left(\frac{M}{10^9
M_{\bigodot}}\right)^{\frac{1}{8}+\frac{3\alpha}{4}}\left(\frac{\chi}{1^{\circ}}\right)^{-1/2}
(k\kappa)^\frac{3}{4}.
$$

The ratio of luminosities in photons, $L_{\gamma }$, and in cosmic rays,
$L_{\rm CR}$, may be estimated by comparing the total available potential
difference in the acceleration region along the rotation axis to its
fraction, spent on the particle acceleration:
$$
\eta=\frac{L_{\gamma}}{L_{\rm CR}}=\frac{{\cal
E}_{\rm max}}{{\cal E}_{\rm curv}},
$$

\begin{equation}
\eta=3.12
\left(\frac{M}{10^9M_{\bigodot}}\right)^{\frac{1}{8}+\frac{3\alpha}{4}}
\left(\frac{\chi}{1^{\circ}}\right)^{1/2}
\frac{Z^{5/4}}{A}
\xi\kappa^{-1/4}k^{3/4},
\label{Eq:eta}
\end{equation}
where $(\xi\kappa^{-1/4})\sim (0.1-2)\ $ for $\ 0<a<M\ $ and $\
R_{\rm S}<H<6R_{\rm S}$. In numerical calculations presented below, we use
$\kappa^{1/4}=0.7$ and $\xi=0.25$, cf. Figs.~\ref{f_1},~\ref{f_2}.

To summarize, the model we use assumes a monochromatic spectrum  of
accelerated particles with $\mathcal{ E}=\mathcal{ E}_{\rm curv}$,
Eq.~(\ref{Eq:Ecurv}), in each particular source. The value of $\mathcal{
E}_{\rm curv}$ depends, within the assumed magnetic-field model, on the
SMBH mass $M$ only (in what follows, we do not consider acceleration of
other particles than protons). The overall flux from the source remains a
free parameter.

\section{Population of the sources and the observed spectrum}
\label{sec:population}
Having discussed the model of particle acceleration in a particular
source, we switch now to the population of sources. As we have seen, the
properties of a single source are determined, within the magnetic-field
model we choose to study, by the SMBH mass (the dependence from the spin
is weak). To reconstruct the UHECR spectrum one has to consider the
population of SMBHs distributed in mass and luminosity. For simplicity we
will assume that the \textit{mean} SMBH luminosity in cosmic rays is
related to its mass,
\begin{equation}
L_{\rm CR} \propto M^{\beta},
\label{Eq:Lcr}
\end{equation}
where $\beta$ is an additional model parameter.
Note that
not every black hole can work as a source, because the
source should possess some special properties (for example, small
inclination angle, the vacuum gap larger than $R_{\rm S}$,
 absence of  numerous charged particles in the vicinity of the black hole
which might imply a thin or even absent
accretion disc). The fraction of sources where the mechanism works is also
encoded in the mean luminosity, Eq.~(\ref{Eq:Lcr}). The observed
spectrum can be obtained by convolving the SMBH mass function with the
(monochromatic) single-source spectrum and Eq.~(\ref{Eq:Lcr}) and taking
into account the propagation effects.

\subsection{The SMBH mass function}
\label{sec:pop:mass-function}
Since dynamical measurements of SMBH masses are available for a very
limited number of cases only, it is a difficult task to find the mass
function precisely. Fortunately, a number of indirect methods to estimate
the SMBH mass are available (see e.g.\ Ref.~\cite{FF} for a review).
Despite having large uncertainties in individual measurements, these
methods are suitable for obtaining average characteristics of the SMBH
population, see e.g.\ Ref.~\cite{MassFunctionReview} for a review.

An important feature of the SMBH mass function is its evolution: SMBHs
grow fast, cf.\ e.g.\ Ref.~\cite{MassFunctionGrow}. For our calculation,
we use one of the most recent published redshift-dependent mass functions
\cite{MassFunctionWeUse}. Of two functions presented there, we choose to
use the one based on the stellar mass functions because it has smaller
statistical uncertainties. The systematic uncertainties of the mass
function may be judged from Ref.~\cite{MassFunctionReview} and are well
within the overall precision of our toy model.

\subsection{The observed spectrum}
\label{sec:pop:spectrum}
Before reaching the Earth, the accelerated protons may interact with
the cosmic microwave background. The main two processes modifying the shape
of the propagated cosmic-ray spectrum are photopion production and $e^+
e^-$-pair production. The former leads to a strong suppression of
the proton flux above few tens of EeV known as the GZK effect~\cite{G,
ZK},
why the latter mostly dominates the attenuation below 10 EeV leading to
the so-called ``dip'' feature in the spectrum~\cite{Hill:1983mk,
Berezinsky:1988wi}. We use the numerical code developed in
Refs.~\cite{propagation}. The code also traces secondary particles
produced
in the interactions. It makes use of the kinetic-equation approach and
calculates the propagation of  nucleons, stable leptons  and photons using
the standard dominant processes (see e.g. Ref.~\cite{reviewsPropag}).
\begin{figure}
\centering
\includegraphics[width=0.7\columnwidth,clip=true,angle=270]{AugerHibridMult1.ps}
\caption{The cosmic-ray flux predicted by the model for three different
assumptions about the SMBH magnetic field, see labels on the plot and
explanations in the text, versus the Auger experimental
data~\cite{AugerSpectrum}.}
\label{f_SpecAuger}
\end{figure}

Fig.~\ref{f_SpecAuger}
presents the predicted cosmic-ray
fluxes in the best-fit model for the
Auger spectrum \cite{AugerSpectrum} for different dependencies of
the magnetic field $B_0$ on SMBH mass mentioned in previous section. The
red curves correspond to $B_0$ given by the Eddington limit, the green one
describes the Shakura--Syunyaev model~\cite{ShakuraSyunyaev,
NovikovThorne}
and the blue one corresponds to the model of Ref.~\cite{OtherBmodel}. The
overall flux normalization is a free parameter. Besides we tried two values
of angle $\chi=1^{\circ}$ and $5^{\circ}$ and varied the
luminosity dependence~(\ref{Eq:Lcr}) parameter $\beta$ in the range
$-1<\beta<2$. The best-fit parameter values are indicated on the plot.
One can see that the first two models produce satisfactory spectral fits
above 10~EeV.

\section{Constraints}
\label{sec:constraints}
In this section, we discuss additional consistency
checks of the model. They include estimates of the accompanying gamma
radiation which should not be in conflict with the measured diffuse
gamma-ray background, estimates of the concentration of sources and of the
luminosity of a single source. We will see that the model passes these
tests. For order-of-magnitude estimates in this section, we assume
$B_{0}\sim B_{\rm Ed}$.

\subsection{Concentration of the sources}
\label{sec:constraints:density}
Let us check that the local concentration
of sources of cosmic rays with energies ${\cal E}\gtrsim  6\times
10^{19}~\mathrm{eV}$  is not in conflict with the lower
limit~\cite{PAO-clustering} based on the statistics of clustering. To this
end, we integrate the SMBH mass function over the range of masses
corresponding to these energies.

The dependence of the particle energy from the black-hole mass is given by
Eq.~(\ref{Eq:EfromMass}). Every black-hole mass corresponds to a range of
particle energies due to the variations in the value of $\kappa$. The
black-hole mass function includes all the black holes with fixed masses
and so all the black holes with every value of $\kappa$. Thus, the lower
limit of the required mass interval is determined by ${\cal E}= 6\times
10^{19}~\mathrm{eV}$ and the maximum value of $\kappa$. It was shown  in
Ref.~\cite{NerSemTk} that particles can be accelerated to the energies of
about $10^{20}~\mathrm{eV}$ only if the size of the vacuum gap is not
smaller than the Schwarzschild radius. In this case, the highest value of
$\kappa^{1/4}$ is $\sim 0.7$ (see Fig.\ref{f_2}).
\begin{figure}
\centering \includegraphics[width=0.48\textwidth]{Energy1.eps}
\caption{\label{f_2}
Dependence of the coefficient $\kappa^{1/4}$ from the angular
momentum $a$ of the black hole with mass $M$ for different values of the
vacuum-gap size $H$ (curves labeled by values of $H/M$; see the text for
definitions of parameters).
}
\end{figure}
Black-hole masses corresponding to particle energies  ${\cal E}\gtrsim
6\times 10^{19}~\mathrm{eV}$ are then $M\gtrsim 10^7 M_{\bigodot}$.
Integrating the mass function  in this range of masses, we obtain
$$
n=\int_{M_{\rm min}}^{M_{\rm max}}\frac{dn}{d\log M}d\log M.
$$
The
integral is saturated at its lower limit that is $M_{\rm max}$ can be
taken arbitrary high to obtain the following estimate,
\begin{equation}
\label{3}
n\sim  10^{-3}\ \frac{1}{\mathrm{Mpc}^3}.
\end{equation}
As we see, the total concentration of sources is larger than the
clustering lower bound of $10^{-4}\frac{1}{\mathrm{Mpc}^3}$. As we have
discussed above, only a
fraction of the calculated concentration $n$ corresponds to the true
concentration of the sources. Our estimate tells us that this fraction
should be not less than a few per cent which is reasonable.

\subsection{Luminosity of a single source}
\label{sec:constraints:luminosity}
A simple estimate of the luminosity of a single source may be obtained as
follows. Consider the observed
flux of cosmic rays with energies ${\cal E}\gtrsim  6\times
10^{19}~\mathrm{eV}$.  The value $j({\cal E})$, which is often reported, is
\begin{equation}
j({\cal E})=
\frac{1}{4\pi} \frac{1}{{\cal
E}}\frac{dF}{d{\cal E}}=\frac{1}{4\pi}\frac{dN}{d{\cal E}},
\label{Eq:j}
\end{equation}
where $N({\cal E})$ is the number of particles with energy $\leq {\cal E}
$, per unit area per unit time. Then, for the flux $F$ we have
$$
{\cal E}\cdot j({\cal E})=\frac{1}{4\pi}\frac{dF}{d{\cal E}},
$$
and, from the recent data~\cite{PAO-GZK, TA-GZK},
we obtain the estimate $F\sim 1.3 \times 10^5 \frac{\rm eV}{\rm m^2\cdot
s}. $

On the other hand, this flux is produced by the sources situated in the
GZK sphere, where the GZK horizon radius is $R_{\rm GZK}\sim 130\
 \mathrm{Mpc}$ for energies ${\cal E}\gtrsim  6\times 10^{19}~\mathrm{eV}$,
 see e.g.~\cite{GZKhorizon}.
We assume that all these sources have
approximately  equal cosmic-ray luminosities $L_0 \left[\frac{\rm eV}{\rm
s}\right]$, independent from the source black-hole mass (in reasonable
agreement with the best-fit values of $\beta$, Sec.~\ref{sec:population}).
Because in this case the GZK radius corresponds to very small redshifts
$z<0.1$, we can neglect changing of the source concentration with the
redshift. For the order-of-magnitude estimate, we neglect also the
difference between the energy with which the particle was emitted and the
final particle's energy with which we detect it on the Earth. Let us also
note that, because we are interested in all sources of the accelerated
 particles with energies ${\cal E}\geq 6\times 10^{19} \mathrm{eV}$, in
the expression for the total flux we have to substitute the total
concentration of sources, Eq.~(\ref{3}), that was obtained by integrating
the mass function.

The flux from every single source, situated at the distance $d$ from us,
is $F_{0}=L_{0}/(4\pi d^2)$. Thus for the total flux we have
$$
F=
\int_0^{R_{\rm GZK}} F_{0}\ n \ dV=
L_{0}\cdot  n \cdot R_{\rm GZK}.
$$
Using
 $
 n\sim 10^{-3}\
 \frac{1}{\mathrm{Mpc}^3}$,
 we obtain
$
L_{0}\sim
6\times 10^{39}\
\frac{\mathrm{erg}}{\mathrm{s}}.
$
The corresponding luminosity in photons  is
$
L_{\gamma}=\eta\cdot L_{0},
$
where $\eta$ is given by Eq.~(\ref{Eq:eta}).
One has
\begin{equation}
\label{4}
L_{\gamma}\sim  10^{40}\
\left(\frac{M}{10^{9}M_{\bigodot}}\right)^{1/8}
{\left(\frac{\chi}{1^{\circ}}\right)}^{1/2}\
\frac{\mathrm{erg}}{\mathrm{s}} .
\end{equation}
This value is much smaller than the typical bolometric luminosity of an
AGN,
$L_{\rm AGN}\sim \left( 10^{41}-10^{43} \right) \
\frac{\mathrm{erg}}{\mathrm{s}}$,
and by far does not exceed the Eddington limit,
 $ L_{\rm Ed}=10^{47} \ \left(\frac{M}{10^{9}M_{\bigodot}}\right)\
\frac{\mathrm{erg}}{\mathrm{s}}$.
Taking the concentration  of the sources of order of the lower limit,
$n\sim 10^{-4}\frac{1}{\mathrm{Mpc}^3}$, does not result in a conflict as
well.

In all the cases, the luminosity is not that far
from the luminosity of an AGN, and a natural question
arises: can we see our sources in $\mathrm{TeV}$ as  point sources? A
simple estimate of the flux of $\mathrm{TeV}$ photons associated with
the particle acceleration from one of the nearest sources, $R\sim
10~\mathrm{Mpc}$, gives the answer: we cannot see them as point sources
because the flux from a single source is smaller than the sensitivity of
the telescopes. Indeed, the flux of $ F_{\gamma}\sim 10^{-1}
\frac{\rm eV}{\mathrm{cm}^2\cdot\mathrm{s}} $ corresponds to the counting
rate of $ \sim 10^{-13}\frac{1}{\mathrm{cm}^2\cdot\mathrm{s}}, $ beyond
the reach of current TeV telescopes. Of course, this does not mean that
strong TeV sources cannot accelerate UHE particles by this mechanism: the
TeV emission may have a totally different origin.

\subsection{Diffuse gamma-ray background}
\label{sec:constraints:diffuse}
While individual cosmic-ray sources have quite low gamma-ray luminosities,
one may wonder about the total emission of all sources in the Universe
(beyond the GZK sphere). The emitted curvature photons have energies of
order a few TeV and interact with the infrared background radiation to
produce electromagnetic cascades in which the energy of the leading gamma
rays downgrade to the GeV band. Electrons in the cascade are deflected by
cosmic magnetic fields so distant sources contribute to the diffuse
gamma-ray background. Let us check that this contribution does not exceed
the measured value of the diffuse flux.

A simple estimate may be obtained as follows.
Consider particles with energies
${\cal E} \sim 10^{20}\mathrm{eV}$.
Cosmic rays with these energies arrive to us
from the interior of the GZK sphere only, but
the associated
photons come from all sources at all distances. Knowing how
much cosmic rays come from the interior of the GZK sphere, we can estimate
how much of them are present in the Universe (keeping in mind that the
number of sources depends on the distance). And supposing that the
luminosity of a source in cosmic rays is connected with its luminosity in
associated photons, Eq.~(\ref{Eq:eta}), we can estimate the total photon
emission from all the sources.

The flux at the Earth, Eq.~(\ref{Eq:j}), is expressed in terms of the
energy at detection, $\cal E$, which in general differs from the energy at
injection, $E_{\rm in}$. In particular, the account of the Universe
expansion (even neglecting additional energy losses) results in ${\cal
E}=E_{\rm in}(1+z)^{-1}$, where
$z$ is the redshift of the source.
The number of the emitted particles per unit time is also $(1+z)$ times
higher than the number of detected particles on the Earth per unit time.

The contribution to $j(E)$ from the sources located at redshift
$z$ is
$$
dj(E,z)=\frac{1}{4\pi}\ \frac{dn_{\rm BH}(M(E_{\rm in}))}{dE_{\rm in}} \
\frac{dN_{0}}{dt}\ \frac{1}{S(z)}\,dV(z),
$$
where $n_{\rm BH}(M(E_{\rm in}))$ is the number density of black holes with
masses $\leq M(E_{\rm in})$ (for $M(E)$, see Eq.~(\ref{Eq:EfromMass})) at
the redshift $z$, $dN_{0}/dt$ is a number of detected particles per unit
time from one of the sources, which were emitted with the energy $E_{\rm
in}$, $dV$ is the volume of a spherical layer at the distance $z$ from us,
$S$ is the area of a sphere with the radius equal to the distance from the
source to the Earth (it is necessary for calculating a number of particles
through the unit area on the Earth). Taking into account the dependence
from the redshift $z$, we may find the source cosmic-ray luminosity
$L(E_{\rm in})$, that is  its total energy emission in cosmic rays per
unit time,
$$
L(E_{\rm in})=\frac{dE_{\rm tot}}{dt}=E_{\rm in}\frac{dN_{0}^{\rm in}}{dt},
$$
where
$dN_{0}^{\rm in}/dt$ is a number of particles with the energy $E_{\rm in}$ emitted
by the source per unit time
and we assumed
that the source emits particles with
the only energy, determined by its mass $M$.

For the flat Universe, we have
$
S=4\pi a_0^2
\rho^{2}(z)$ and $ dV=Sa_{0}d\rho(z),
$
where
$$
 \rho(z)=\int_0^z\frac{dz}{H(z)a_{0}}
 $$
 is the geodesic coordinate distance
 from the
observer to the source location.
Finally
\begin{widetext}
 \begin{equation}
\label{f_9}
j(E)=
\int_0^{z_{\rm max}} \!
\frac{1}{4\pi}\frac{dn_{\rm BH}(M((1+z)E))}{d\log
\frac{M}{M_{\bigodot}}} \frac{L((1+z)E)}{(1+z)^{2}E}\ \left.
\frac{d\log \frac{M}{M_{\bigodot}}}{dE_{\rm in}}\right|_{M((1+z)E))}
 a_{0}d\rho(z),
\end{equation}
\end{widetext}
where $dn_{\rm BH}(M)/d\log M$ is just the mass function for a given
$z$~\cite{MassFunctionWeUse},
$z_{\rm max}$ is the redshift of the most
distant source.
Here we are interested in the values of the mass function at the points
$M((1+z)E)$ as a function of z. This function could be easily constructed
using the data from Ref.~\cite{MassFunctionWeUse}. Using
Eq.~(\ref{Eq:EfromMass}), we have
$$
\frac{d\log
\frac{M}{M_{\bigodot}}}{dE_{\rm in}}(M((1+z)E))
\sim \frac{1.16}{E_{\rm in}}=\frac{1.16}{(1+z)E}.
$$

For simplicity, let us suppose  that all the sources have the same
luminosities $L((1+z)E)=L_{0}$. This assumption is in reasonable agreement
with the results obtained in section~\ref{sec:population}. Thus we can
take the luminosity out of the integral. Let us now imagine just for a
moment, that cosmic rays with such energies could come from all distances,
and calculate the ratio of the values $j(E)$ for the cosmic rays with the
energies $\sim 10^{20}$~eV from the GZK sphere ($z_{\rm max}=z_{\rm
GZK}\sim 0.01$) and from the whole Universe (e.g.\ $z_{\rm max}=2$). The
luminosities in front of these two integrals are cancelled and after
integrating we obtain:
$$
\frac{E^2 j(E)_{\rm tot}}{E^2
j(E)_{\rm GZK}}=20.
$$
Taking into account the observed value of the cosmic-ray flux, we estimate
$$
E^2 j(E)^{\rm tot}\sim 10^5 \frac{\mathrm{eV}}{\mathrm{m}^2
\cdot\mathrm{s}\cdot \mathrm{sr}}.
$$

For calculating the value of $E_{\gamma}^2 j(E)_{\gamma}$ for the gamma
radiation from all the possible sources in the whole Universe we should
only replace $d\log M/dE$ in the integral (\ref{f_9}) by the
corresponding expression for photons, Eq.~(\ref{4}),
and the cosmic-ray luminosity $L$ by the photon luminosity $L_{\gamma}=\eta
L$, Eq.~(\ref{Eq:eta}). The order-of-magnitude estimate then reads
$$
E_{\gamma}^2
j(E)_{\gamma}^{\rm tot}=\eta E^2 j(E)^{\rm tot}\sim
10^5\frac{\mathrm{eV}}{\mathrm{m}^2 \cdot\mathrm{s}\cdot \mathrm{sr}}.
$$
This total gamma-ray emission associated
with particle acceleration in all the sources does not exceed the observed
value of the diffuse gamma-ray background \cite{FERMIdiffuse}.

We have also performed a more detailed numerical simulation of the
secondary gamma-ray flux. The injection spectrum of the curvature photons
is similar to the synchrotron one ~\cite{Electrod},

$$
I(\omega)\approx 2\sqrt{3} Z^2 e^2 \gamma
\frac{\omega}{\omega_{c}}\int_{2\omega/\omega_c}^{\infty}K_{5/3}(x)dx,
$$
where
$$
\omega_{c}= 3 \frac{\gamma^3}{R}
$$
is the ``critical'' frequency
(for higher frequencies, radiation is negligible), $\omega$ is the
frequency of the radiated photons, $K_{5/3}$ is the Macdonald function,
$\gamma=({\cal E}/m)$ is the particle Lorentz factor, $R$ is the curvature
radius of the particle trajectory (for curvature radiation, $R$ is
constant and it is equal to the curvature radius of the magnetic lines).
The total flux of the photons from a single source is related to the
cosmic-ray flux by the coefficient $\eta$, Eq.~(\ref{Eq:eta}). For the
best-fit spectrum of Sec.~\ref{sec:pop:spectrum}, we used the same code to
describe the propagation of the accompanying gamma rays and to calculate
the observed gamma-ray flux. The result is presented in
Fig.~\ref{f_SpecAugerGamma}.
\begin{figure}
\centering
\includegraphics[width=0.7\columnwidth,clip=true,angle=270]{AugerHibridMult1all.ps}
\caption{Gamma-ray fluxes predicted by the same models as shown in
Fig.~\ref{f_SpecAuger}.}
\label{f_SpecAugerGamma}
\end{figure}
One may see that the diffuse gamma-ray upper limit is satisfied.

\section{Conclusions and discussion}
\label{sec:concl}
We have constructed and studied a toy model of UHECR acceleration in the
vicinity of numerous and various supermassive black holes in centers of
galaxies. The model assumes that:
\begin{itemize}
 \item
cosmic-ray particles are accelerated by the regular electric field within
a few $R_{S}$ from the SMBH~\cite{NerSemTk}; the field configuration is
given by the solutions of Refs.~\cite{BHsolution1, BHsolution2} and is
fully determined by the SMBH mass $M$, its angular momentum $a$ and the
magnetic-field normalization $B_{0}$;
\item
all cosmic-ray particles accelerated near a given SMBH have similar
initial conditions and therefore all are accelerated up to one and the
same energy limited by the curvature-radiation losses; this maximal energy
is calculated in the model and depends on $M$ only (provided $B_{0}$ is a
given, model-dependent, function of $M$; dependence from $a$ smooths this
monochromatic spectrum insignificantly);
\item
the \textit{mean} flux of a source (which accounts for the fraction of the
sources where this mechanism does work) depends from $M$ in a power-like
manner; the normalization and the exponent are two free parameters of the
model (the best fit to the cosmic-ray spectrum indicates that this
dependence is weak);
\item
the concentration of sources is determined by the redshift-dependent SMBH
mass function taken from astrophysical literature.
\end{itemize}
Within these assumptions and given the $B_{0}(M)$ relation is fixed (we
considered three popular choices for it), the model has two free
parameters which we find by fitting the cosmic-ray spectrum at the Earth
to the experimental data. With parameters fixed in this way, we subject
the model to several further tests which it passes succesfully:
\begin{enumerate}
 \item
the concentration of sources is large enough to satisfy the constraints
from absence of clustering in UHECR arrival directions;
\item
the luminosity of a particular source, determined by the flux
normalization and concentration, is not too high;
\item
secondary gamma rays from distant sources do not overshy the measured GeV
diffuse gamma-ray background.
\end{enumerate}
Given the success of the toy model, it is interesting to discuss its
possible refinements.
The assumptions we have made within the model are quite robust and
realistic. One subtle point is related to the value of the SMBH angular
momentum $a$ which may vary from one black hole to another. However, these
variations are probably modest given the scaling relation between the mass
and the angular momentum of cosmic black holes proposed in
Ref.~\cite{Nature-momentum}. The precision of predictions may be improved
with more realistic modelling of the acceleration mechanism, in
particular, with account of the charge concentration at the SMBH and in
the acceleration region, of the finite thickness of the accretion disk
etc. Ultimately, this approach might give answer to the question, which
particular SMBHs are strong sources and which are not, thus determining
the (presently free) parameter theoretically. However, this is a
complicated task and is far beyond the scope of the present work.

One possible question concerns the low-energy part of the spectrum where,
as is clearly seen from Fig.~\ref{f_SpecAugerGamma}, the contribution of
the mechanism we discuss is insufficient to explain the observed spectrum
due to the depletion of the SMBH mass function at low masses. It is
tempting to speculate that this depletion is compensated by a huge
contribution of the SMBH in our own Galaxy which, indeed, has the
appropriate mass. A quantitative analysis of this proposal requires,
however, a much more precise study of physical properties and
possibilities for particle acceleration close to the Galactic Center.

{\bf Acknowledgements. } We are indebted to V.~Alba, P.~Dunin-Barkowski,
G.~Farrar, A.~Neronov, D.~Semikoz and I.~Tkachev for interesting
discussions. We thank the authors of Ref.~\cite{MassFunctionWeUse} for
providing numerical values of the SMBH mass function we used. This work
was supported in part by the RFBR grant 10-02-01406 and by the grant of
the President of the Russian Federation NS-5590.2012.2, by the RFBR grants
12-02-01203 and 11-02-01528 (S.T.) and by the Dynasty Foundation (K.P.\
and S.T.). Numerical part of the work was performed at the cluster of the
Theoretical Division of INR RAS.