\section{Introduction}

In chemical kinetics, the temporal evolution of a well-stirred system of chemically reacting molecules is classically described using reaction rate equations. Reaction-rate equations are a mean-field description formulated as coupled ordinary differential equations. The number of molecules is continuous in time, and reaction rates are quantified using macroscopic rate constants. Such reaction rate equations, however, do not always provide an accurate description. This is the case especially, but not only, when the number of molecules of the various chemical species (henceforth called the population) is much smaller than Avogadro's number~\cite{Qian:2002,Shibata:2004}. At low population, the number of molecules is not large enough for fluctuations to be negligible. In addition, fluctuations may play an important role in the kinetics~\cite{Qian:2002,Li:2008a}. Even at high population, correlated fluctuations can cause the mean to behave in a way that is not captured by a mean-field description~\cite{Shibata:2004,Gardiner:1976,Engblom:2006}. These effects can be accounted for by stochastic kinetic models, which can incorporate thermal fluctuations in a number of ways. An approach that has become canonical is the chemical master equation (CME)~\cite{Gillespie:1976,Gillespie:1977,Gillespie:1992}, a Markov-chain model with many applications in physics, chemistry, and biology. The CME models the kinetics of any chemical reaction system that is well stirred and thermally equilibrated \cite{Gillespie:1992}. Its high dimensionality, however, renders analytical approaches intractable. \\

Numerical methods to sample trajectories from the CME mostly rely on kinetic Monte Carlo approaches~\cite{Bortz:1975}. The canonical kinetic Monte Carlo approach for sampling a trajectory of the CME is Gillespie's stochastic simulation algorithm (SSA)~\cite{Gillespie:1976,Gillespie:1977,Gillespie:1992}. SSA is governed by the joint probability density 
\begin{equation}
p(\tau,\mu | \vec{n}(t)) = (a \text{e}^{-a\tau})( a_{\mu}/a)
\label{eq:jointpdf}
\end{equation}
for the random variables $\tau$ (the time to the next reaction) and $\mu$ (the index of the next reaction). The vector \mbox{$\vec{n}(t)=(n_{1},\dots,n_{N})$} is the population at time $t$. Each entry \mbox{$n_i$} is the number of molecules of the respective species S\mbox{$_i$}, and $N$ is the total number of species. The propensity of each reaction $\mu$ is defined as \mbox{$a_\mu = c_\mu h_\mu$}, where \mbox{$c_{\mu}$} is the specific probability rate, and \mbox{$h_{\mu}= h_{\mu}(\vec{n})$} is the reaction degeneracy, which is the number of possible combinations of reactants in reaction $\mu$ given the population $\vec{n}$. 
The reaction propensity is such that \mbox{$a_{\mu}\md t$} is the probability that a randomly selected combination of reactant molecules of reaction $\mu$ at time $t$ will react in the next infinitesimal time interval \mbox{$\md t$}. The total propensity is \mbox{$a=\sum_{\mu=1}^{M}a_{\mu}$}, where $M$ is the total number of reactions. \\

Existing SSA formulations can be classified into {\em exact} and {\em approximate} methods. Exact methods sample from the probability density in Eq.~\ref{eq:jointpdf}. These formulations include the direct method (DM)~\cite{Gillespie:1976,Gillespie:1992}, the first reaction method (FRM)~\cite{Gillespie:1976}, Gibson-Bruck's next-reaction method (NRM)~\cite{Gibson:2000}, a Gibson-Bruck variant of the DM~\cite{Gibson:2000}, the optimized direct method (ODM)~\cite{Cao:2004}, the sorting direct method (SDM)~\cite{McCollum:2006}, the logarithmic direct method (LDM) (unpublished, \cite{Li:2006b}), and the composition-rejection formulation (SSA-CR) \cite{Slepoy:2008}. 
Approximate SSA formulations provide better computational efficiency for large numbers of molecules by sampling from an approximation to the probability density in Eq.~\ref{eq:jointpdf}. These methods include $\tau$-leaping~\cite{Gillespie:2001,Cao:2005,Cao:2006,Peng:2007a}, $k_\alpha$-leaping~\cite{Gillespie:2001}, $R$-leaping~\cite{Auger:2006a}, $L$-leap~\cite{Peng:2007}, $K$-leap~\cite{Cai:2007}, the slow-scale method~\cite{Cao:2005a}, and implicit $\tau$-leaping~\cite{Rathinam:2003}.\\ 

In this paper we focus on exact methods. They offer the advantage of being parameter-free, whereas all approximate methods contain parameters that need to be adjusted by the user. The computational cost of exact SSA formulations is dominated by the steps needed to sample the next reaction and update the propensities after a reaction has fired \cite{Gibson:2000,Cao:2004,Slepoy:2008}. In Gillespie's original DM and FRM, this leads to a computational cost that scales linearly with the number of reactions in the network. Various improved SSA formulations have been proposed in order to reduce this computational cost. 
The most notable improvements include the use of dependency graphs to reduce the number of propensities that need to be updated \cite{Gibson:2000}, and various sampling schemes of higher efficiency \cite{Gibson:2000,Cao:2004,McCollum:2006,Slepoy:2008}. All of these sampling schemes can be interpreted as instances of the random-variate generation problem \cite{Slepoy:2008} as described in Devroye's compendium \cite{Devroye:1986} and can reduce the computational cost (CPU time) of sampling the next reaction. These improvements have reduced the computational cost of SSA to logarithmic or even constant scaling for weakly coupled networks. 
For strongly coupled networks, however, the computational cost of all improved SSA formulations still scales linearly with the number of reactions. We define {\em weakly coupled} networks as those where the maximum number of reactions that are influenced by any other reaction, i.e.~the maximum degree of coupling of the network, is independent of system size. This is in contrast to {\em strongly coupled} networks, where the number of influenced reactions grows proportionally with system size and can be as large as the total number of reactions. In such networks, the total number of reactions grows faster than the number of species when the latter is increased. Strongly coupled networks frequently occur, e.g., in nucleation-and-growth models, scale-free biochemical networks, and colloidal aggregation systems. In these cases, the scaling of the computational cost of most of the improved SSA formulations with system size is equivalent to that of DM (see Sec.~\ref{sec:oldssa}). \\

We present a novel SSA formulation with a computational cost that scales at most linearly with the number of {\em species}, making it especially efficient for strongly coupled networks. This is made possible by restricting the class of systems to networks containing only elementary chemical reactions, where every reaction has at most two reactants \cite{Gillespie:1992}. This allows factoring out one of the species from every reaction propensity, leading to {\em partial propensities} that depend on the population of at most one species. Any non-elementary reaction can always be broken down into elementary reactions, at the expense of an increase in system size \cite{Gillespie:1992,Wilhelm:2000,Schneider:2000}. The use of partial propensities leads to SSA formulations with a computational cost that scales as some function of the number of species rather than the number of reactions. \\

In Sec.~\ref{sec:methods_description}, we formally introduce the concept of partial propensities and present two partial propensity variants of the exact SSA: the partial-propensity direct method (PDM) and the sorting partial-propensity direct method (SPDM). They use partial propensities and efficient data structures for sampling the next reaction and for updating the partial propensities after a reaction. We benchmark them in Sec.~\ref{sec:benchmarks} and show that their computational cost scales at most linearly with the number of species in the network, irrespective of the degree of coupling. The benchmarks include two strongly coupled networks, for which the degree of coupling grows with system size, a weakly coupled reaction network with a constant maximum degree of coupling, and a small, fixed-size biological multiscale (stiff) network. In order to test the competitiveness of our algorithm in cases where several other SSA formulations might be more efficient, we choose the most weakly coupled network possible, the linear chain model, where the number of reactions scales linearly (with a proportionality constant of 1) with the number of species \cite{Albert:2005,Cao:2004}. The multiscale biological network is included in order to benchmark the new algorithms on small systems and when the reaction propensities span several orders of magnitude. In Sec.~\ref{sec:conclusion} we summarize the main results, discuss the limitations of the presented method, and give an outlook on possible future developments and applications.

\section{Computational cost of previous exact SSA formulations}\label{sec:oldssa}

We review the scaling of the computational cost of previous exact SSA formulations. In order to express scaling with system size $x$, we use the Bachmann-Landau notation, writing $C(x)\in O(f(x))$ ($C(x)$ is $O(f(x))$) whenever $C(x)>0$ is bounded from above by $f(x)$ as $C(x)\le \alpha f(x)$, for all $x$ and some constant pre-factor \mbox{$\alpha>0$} that is independent of $x$. \\

Since DM and FRM form the basis for most exact SSA's, we first focus on these two. DM's computational cost is \mbox{$O(M)$}~\cite{Gillespie:1976,Gillespie:1992,Gibson:2000,Cao:2004,Slepoy:2008}, where $M$ is the total number of reactions (see also Appendix \ref{app:ssa}). In FRM, the sampling strategy for $\mu$ and $\tau$ is different (see Appendix \ref{app:ssa}). This, however, does not change the scaling of the computational cost of FRM, which remains \mbox{$O(M)$}. Since the FRM sampling strategy involves discarding \mbox{$M-1$} reaction times, its computational cost generally has a larger pre-factor than that of DM~\cite{Gillespie:1976,Gibson:2000,Cao:2004}. \\

NRM is an improvement over FRM in which the \mbox{$M-1$} unused reaction times are suitably reused, and data structures such as indexed priority queues and dependency graphs are introduced. The indexed priority queue, which is equivalent to a heap tree, is used to sort the \mbox{$\tau_i$}'s more efficiently; the dependency graph is a data structure that contains the indices of the reactions whose propensities are to be recomputed after a certain reaction $\mu$ has fired. This avoids having to recompute all \mbox{$a_\mu$}'s after every reaction. Each reaction is represented as a node in the dependency graph, and nodes $i$ and $j$ are connected by a directed edge if and only if the execution of reaction $i$ affects the propensity (through the population of reactants) of reaction $j$. These data structures, together with the reuse of reaction times, reduce the computational cost of NRM to \mbox{$O(k\log _2 M)$}, where $k$ is the out-degree of the dependency graph, that is, the degree of coupling of the reaction network. In strongly coupled networks, $k$ is a function of $M$ and is $O(M)$. The computational cost of NRM is thus \mbox{$O(M)$} for strongly coupled networks. Even for some weakly coupled networks, the computational cost of NRM has been empirically shown to be \mbox{$O(M)$}~\cite{Cao:2004}. This is due to the additional overhead, memory-access operations, and cache misses introduced by the complex data structures (indexed priority queue, dependency graph) of NRM. The scaling of the computational cost of the Gibson-Bruck variant of DM is equal to that of NRM, albeit with a larger pre-factor~\cite{Gibson:2000}. For weakly coupled networks where $k(M)$ is $O(1)$, independent of system size, the computational cost is further reduced to $O(1)$ in the SSA-CR formulation \cite{Slepoy:2008} under the assumption that the ratio of maximum to minimum propensity is bounded. For strongly coupled networks, where $k(M)$ is $O(M)$, the computational cost of SSA-CR is $O(M)$ \cite{Slepoy:2008,Schulze:2008}. \\

ODM is an improvement over DM where the reactions are sorted in descending order of firing frequency. This makes it more probable to find the next reaction close to the beginning of the list and, hence, reduces the search depth for finding the index of the next reaction using linear search. ODM estimates the firing frequencies of all reactions during a short pre-simulation run of about 5--10\% of the length of the entire simulation \cite{Cao:2004,McCollum:2006}. In order to reduce the cost of updating the propensities after a reaction has fired, ODM also uses a dependency graph. Irrespective of the degree of coupling, the computational cost of ODM is \mbox{$O(M)$}, which was confirmed in benchmarks \cite{Cao:2004}. SDM is a variant of ODM that does not use pre-simulation runs, but dynamically shifts up a reaction in the reaction list whenever it fires (``bubbling up'' the more frequent reactions). This further reduces the pre-factor of the computational cost of SDM compared to that of ODM, but the scaling remains \mbox{$O(M)$} \cite{McCollum:2006}. \\

LDM uses a binary search tree (recursive bisection) on an ordered linear list of cumulative sums of propensities to find the next reaction. This is reported to reduce the average search depth of this step to \mbox{$O(\log _2 M)$}~\cite{Li:2006b}. Irrespective of the degree of coupling, however, the update step is $O(M)$ since on average \mbox{$(M+1)/2$} sums of propensities need to be recomputed, rendering the computational cost of LDM \mbox{$O(M)$}. \\

In summary, the computational cost of previously reported exact SSA formulations is \mbox{$O(M)$} for strongly coupled networks. For weakly coupled networks, however, some are significantly more efficient and can be $O(\log _2 M)$ or even $O(1)$. 

\section{Partial-propensity methods}\label{sec:methods_description}

We introduce the concept of partial propensities for elementary reactions and use it to formulate two partial-propensity direct methods, PDM and SPDM, whose computational cost scales at most linearly with the number of species, even for strongly coupled networks. SPDM uses concepts from SDM \cite{McCollum:2006} to dynamically rearrange reactions, which reduces the average search depth for sampling the next reaction in a multiscale network. \\

We define the partial propensity of a reaction \rev{with respect to one of its reactants} as the propensity per molecule of \rev{this reactant}. For example, the partial propensity \mbox{$\pi_{\mu}^{(i)}$} of reaction \mbox{$\mu$} with respect to (perhaps the only) reactant S\mbox{$_i$} is \mbox{$a_{\mu}/n_{i}$}, where \mbox{$a_{\mu}$} is the propensity of reaction \mbox{$\mu$} and \mbox{$n_i$} is the number of molecules of S\mbox{$_i$}. The partial propensities of the three elementary reaction types are:
\begin{itemize}
\item Bimolecular reactions (\mbox{$\r{S}_i\,+\,\r{S}_j\,\rightarrow$} Products): \mbox{$a_{\mu}\,=\,n_i\,n_j\,c_{\mu}$} and \mbox{$\pi_{\mu}^{(i)}\,=\,n_j\,c_{\mu}$}, \mbox{$\pi_{\mu}^{(j)}\,=\,n_i\,c_{\mu}$}. \\ 
\rev{If both reactants are of the same species, i.e.~\mbox{$\r{S}_i = \r{S}_j$}, only one partial propensity exists, $\pi_{\mu}^{(i)} \, = \, \frac{1}{2} (n_i - 1)c_{\mu}$ because the reaction degeneracy is \mbox{$\frac{1}{2} n_i (n_i - 1)$}. If $n_i = 0$, the partial propensity becomes negative. As explained in the caption of Fig.~\ref{Example_PDM} this, however, does not require any special treatment.}
\item Unimolecular reactions (\mbox{$\r{S}_i\,\rightarrow$} Products): \mbox{$a_{\mu}\,=\,n_i\,c_{\mu}$} and \mbox{$\pi_{\mu}^{(i)}\,=\,c_{\mu}$}.
\item Source reactions (\mbox{$\emptyset\,\rightarrow$} Products): \mbox{$a_{\mu}\,=\,c_{\mu}$} and \mbox{$\pi_{\mu}^{(0)}\,=\,c_{\mu}$}.
\end{itemize}    
We consider only these elementary reaction types since any reaction with three of more reactants can be treated by decomposing it into a combination of elementary reactions~\cite{Gillespie:1992,Wilhelm:2000,Schneider:2000}. \\   

\subsection{The partial-propensity direct method (PDM)}\label{sec:PDM_description}

In PDM, the index of the next reaction $\mu$ is sampled in a way that is algebraically equivalent to that of DM, as shown in Appendix \ref{pdm:proof}.
The major novelties in PDM are the use of partial propensities and efficient data structures that reduce the number of operations needed to sample $\mu$ and to update the partial propensities. The time to the next reaction is sampled as in DM. We first present the main principles behind the new sampling and update schemes and then describe them in detail. The complete algorithm is given in Table \ref{our_algorithm}.\\

\subsubsection{Main principles behind PDM}

PDM uses partial propensities and groups them in order to efficiently sample the index of the next reaction and update the partial propensities after a reaction has fired. For the sampling step, the partial propensities are grouped according to the index of the factored-out reactant, yielding at most \mbox{$N+1$} groups of size $O(N)$. Sampling then proceeds in two steps: we first sample the index of the group before sampling the actual partial propensity inside that group. This grouping scheme reduces the number of operations needed for sampling the next reaction using a concept that is reminiscent of two-dimensional cell lists \cite{Hockney:1988}. If all partial propensities are in the same group, or if every group contains only a single partial propensity, the sampling step of PDM is no more efficient than that of DM. These cases, however, can only occur if the function \mbox{$M(N)$} is \mbox{$O(N)$} (see for example the linear chain model) and both PDM and DM hence have a computational cost of $O(N)$ for sampling the index of the next reaction. \\

After the selected reaction has been executed, we use a dependency graph over species (partial propensities), rather than reactions, to find all partial propensities that need to be updated. This is possible because partial propensities depend on the population of at most one species, and is analogous to a Verlet list \cite{Verlet:1967}. This limits the number of updates to be $O(N)$. 
In addition, partial propensities of unimolecular reactions are constant and never need to be updated. 
In weakly coupled networks, where the degree of coupling is $O(1)$, the scaling of the computational cost of the update becomes equal to that of methods that use dependency graphs over reactions, such as SSA-CR, ODM, and SDM. \\

We illustrate the sampling scheme of PDM in a simple protein aggregation example. Consider proteins that aggregate to form at most tetrameric complexes. There are \mbox{$N=4$} species in the reaction network: monomers, dimers, trimers, and tetramers. All species except tetramers can aggregate in all possible combinations to form multimeric complexes (4 bimolecular reactions). In addition, all multimeric complexes can dissociate into any possible combination of two smaller units (4 unimolecular reactions) and monomers are constantly produced (1 source reaction). This reaction network is described by $M=9$ partial propensities $(\pi_{1}^{(0)})$, $(\pi_{2}^{(1)}$, $\pi_{3}^{(1)}$, $\pi_{4}^{(1)})$, $(\pi_{5}^{(2)}$, $\pi_{6}^{(2)})$, $(\pi_{7}^{(3)})$, $(\pi_{8}^{(4)}$, $\pi_{9}^{(4)})$. Grouping the partial propensities according to the index of the factored-out reactant given in the superscript, we obtain 5 (\mbox{$=N+1$}) groups as indicated by the parentheses. Along with each group, we store the sum of all partial propensities inside it. Using a random number, we sample the group that contains the next reaction, before finding the corresponding partial propensity inside that group. Assume that in our example reaction 7 is to fire next. The search depth to find the group index is 4 and we need 1 additional operation to find the partial propensity $(\pi_{7}^{(3)})$. PDM thus requires 5 operations to sample the next reaction in this network of 9 reactions. The average search depth of sampling the next reaction in this example is \mbox{$37/9\approx 4.1$}.  \\      

In the next section, we formally describe PDM and its data structures.     

\subsubsection{Detailed description of the PDM algorithm}

All partial propensities are stored in the ``partial-propensity structure'' \mbox{${\boldsymbol\Pi}=\left\{ \boldsymbol{\Pi}_i \right\}_{i=0}^N$} as a one-dimensional array of one-dimensional arrays \mbox{${\boldsymbol\Pi}_i$}. Each array \mbox{${\boldsymbol\Pi}_i$} contains the partial propensities belonging to group $i$. 
The partial propensities of source reactions are stored as consecutive entries of the 0\mbox{$^\text{th}$} array \mbox{$\boldsymbol \Pi _0$}. The partial propensities of all reactions that have species S\mbox{$_1$} as one of its reactants are stored as consecutive entries of \mbox{$\boldsymbol \Pi _1$}. In general, the \mbox{$i^\text{th}$} array \mbox{$\boldsymbol \Pi_{i}$} contains the partial propensities of all reactions that have S\mbox{$_i$} as a reactant, provided these reactions have not yet been included in any of the previous \mbox{$\boldsymbol \Pi_{j<i}$}. That is, out of the two partial propensities of a reaction \mbox{$\mu$} with S\mbox{$_i$} and S\mbox{$_j$} as its reactants, \mbox{$\pi_{\mu}^{(i)}$} is part of \mbox{$\boldsymbol \Pi_{i}$} if \mbox{$i<j$}, and \mbox{$\pi_{\mu}^{(j)}$} is not stored anywhere. Notice that, since the different \mbox{$\boldsymbol \Pi_{i}$}'s can be of different length, storing them as an array of arrays is more (memory) efficient than using a matrix (i.e.~a two-dimensional array). The reaction indices of the partial propensities in ${\boldsymbol \Pi}$ are stored in a look-up table \mbox{$\textbf{L}=\{\mathrm{\mathbf{L}}_i\}_{i=0}^{N}$}, which is also an array of arrays. This makes every reaction \rev{\mbox{$\mu$}} identifiable by a unique pair of indices, \rev{a group index} $I$ and \rev{an element index} $J$, such that the partial propensity of reaction \mbox{$\mu =\r{L}_{I,J}$} is stored in \mbox{$\Pi_{I,J}$}. \\

We further define the ``group-sum array'' $\boldsymbol \Lambda$, storing the sums of the partial propensities in each group \mbox{$\boldsymbol \Pi_{i}$}, thus \mbox{$\Lambda _i = \sum_{j} \Pi_{i,j}$}, \mbox{$i\,=\,0,\ldots, N$}. In addition, we also define \mbox{$\boldsymbol{\Sigma}$}, the array of the total propensities of all groups, as \mbox{$\Sigma_i = n_{i}\Lambda _i$}, \mbox{$i\,=\,1,\ldots, N$}, and \mbox{$\Sigma _0 = \Lambda _0$}. The total propensity of all reactions is then \mbox{$a=\sum_{i=0}^{N}\Sigma_i$}. The use of \mbox{$\boldsymbol\Lambda$} avoids having to recompute the sum of all partial propensities in \mbox{$\boldsymbol \Pi_{i}$} after one of them has changed. Rather, the same change is also applied to \mbox{$\Lambda _i$} and computing the new \mbox{$\Sigma _i$} only requires a single multiplication by \mbox{$n_i$}. Using these data structures and a single uniformly distributed random number \rev{\mbox{$r_1\,\in\,[0,1)$}}, the next reaction $\mu$ can efficiently be sampled in two steps: (1) sampling the \rev{group} index $I$ such that

\begin{equation}
\label{row_mu}
I ={\operatorname{min}}\left[ I' \,\, : \,\, r_{1}a < \sum_{i=0}^{I'} \Sigma_{i} \right]
\,
\end{equation}
and (2) sampling \rev{the} \rev{element} index $J$ in ${\boldsymbol \Pi}_I$ such that
\begin{equation}
J ={\operatorname{min}}\left[ 
J ' \,\, : \,\, r_{1}a < \sum_{j=1}^{J'}n_{I}\Pi_{I,j} + \left( \sum_{i=0}^{I} \Sigma_{i}\right) - \Sigma_{I}  \right] \,.
\label{col_mu}
\end{equation}
(See Appendix~\ref{pdm:proof} for a proof of the equivalence of this sampling scheme to that of DM.) Using the temporary variables
\begin{equation}
\label{col_mu1}
\Phi = \sum_{i=0}^{I} \Sigma_{i},\qquad
\Psi = \frac{r_1 a - \Phi + \Sigma _I}{n_I} \, ,
\end{equation}
Eq.~\ref{col_mu} can be efficiently implemented as
\begin{equation}
\label{col_mu3}
J ={\operatorname{min}}\left[ 
J' \,\, : \,\, \Psi < \sum_{j=1}^{J'}\Pi_{I,j} \right] \, .
\end{equation}

The indices $I$ and $J$ are then translated back to the reaction index $\mu$ using the look-up table $\mathbf{L}$, thus \mbox{$\mu = \r{L}_{I,J}$}. \\

Once a reaction has been executed, $\vec{n}$, $\boldsymbol \Pi$, $\boldsymbol \Lambda$, and $\boldsymbol \Sigma$ need to be updated. This is efficiently done using three update structures:
\begin{enumerate}
\item[\mbox{$\mathrm{\mathbf{U}^{(1)}}$}] is a array of $M$ arrays, where the \mbox{$i^\text{th}$} array contains the indices of all species involved in the \mbox{$i^\text{th}$} reaction. 
\item[\mbox{$\mathrm{\mathbf{U}^{(2)}}$}] is a array of $M$ arrays containing the corresponding stoichiometry (the change in population of each species upon reaction) of the species stored in \mbox{$\mathrm{\mathbf{U}^{(1)}}$}.
\item[\mbox{$\mathrm{\mathbf{U}^{(3)}}$}] is a array of $N$ arrays, where the \mbox{$i^\text{th}$} array contains the indices of all entries in $\boldsymbol \Pi$ that depend on \mbox{$n_{i}$}, thus: 
\end{enumerate}
\begin{eqnarray}
\mathbf{U}^{(3)}\,=\,
\left\{\begin{array}{l}
\mathbf{U}_1^{(3)} \,=\,\left(i^{1}_1,j^{1}_1 \quad i^{1}_2,j^{1}_2 \quad\ldots\quad\ldots\quad\ldots\,\right)\\
\mathbf{U}_2^{(3)} \,=\,\left(i^{2}_1,j^{2}_1 \quad i^{2}_2,j^{2}_2 \quad\ldots\,\right)\\
\vdots\\
\mathbf{U}_N^{(3)} \,=\,\left(i^{N}_1,j^{N}_1 \quad i^{N}_2,j^{N}_2 \quad\ldots\quad\ldots\,\right).
\end{array}\right.
\label{U3}
\end{eqnarray}

When a reaction is executed, the populations of the species involved in this reaction change. Hence, all entries in $\boldsymbol \Pi$ that depend on these populations need to be updated. After each reaction, we use \mbox{$\mathrm{\mathbf{U}}^{(1)}$} to determine the indices of all species involved in this reaction. The stoichiometry is then looked up in \mbox{$\mathrm{\mathbf{U}}^{(2)}$} and the population $\vec{n}$ is updated. Subsequently, \mbox{$\mathrm{\mathbf{U}^{(3)}}$} is used to locate the affected entries in $\boldsymbol \Pi$ and recompute them. The two data structures \mbox{$\mathrm{\mathbf{U}}^{(1)}$} and \mbox{$\mathrm{\mathbf{U}}^{(2)}$} are a sparse representation of the stoichiometry matrix, and \mbox{$\mathrm{\mathbf{U}}^{(3)}$} represents the dependency graph over species. Since the partial propensities of unimolecular and source reactions are constant and need never be updated, \mbox{$\mathrm{\mathbf{U}^{(3)}}$} only contains the indices of the partial propensities of bimolecular reactions.
The size of \mbox{$\mathbf{U}^{(3)}$} is at most a factor of $N$ smaller than that of the corresponding dependency graph over reactions, since partial propensities depend on the population of at most one species.  
Figure \ref{Example_PDM} summarizes the data structures used in PDM for an example reaction network. The complete algorithm is given in Table~\ref{our_algorithm}. 
Overall, PDM's computational cost is \mbox{$O(N)$} and its memory requirement is \mbox{$O(M)$}, irrespective of the degree of coupling (see Appendix \ref{app:complexity}).  

\subsection{The sorting partial-propensity direct method (SPDM)}

The sorting partial-propensity direct method (SPDM) is the partial-propensity variant of SDM~\cite{McCollum:2006}. In SPDM, the group and element indices $I$ and $J$ are bubbled up whenever the reaction \mbox{$\mu = \r{L}_{I,J}$} fires. The reordered indices are stored in an array for $I$, and an array of arrays of the size of $\boldsymbol \Pi$ for the $J$'s. This requires an additional \mbox{$N+M$} memory, but further reduces the search depth to sample the next reaction, especially in a multiscale (stiff) network. The computational cost of SPDM is also $O(N)$ (see Appendix \ref{app:complexity}), but with a possibly reduced pre-factor. 

\section{Benchmarks}\label{sec:benchmarks}

We benchmark the computational performance of PDM and SPDM using four chemical reaction networks that are prototypical of: (a) strongly coupled reaction networks, (b) strongly coupled reaction networks comprising only bimolecular reactions, (c) weakly coupled reaction networks, and (d) multiscale biological networks. The first two benchmarks consider strongly coupled networks where the degree of coupling scales with system size (see column ``degree of coupling'' in Table \ref{table_test_cases}). The first benchmark consists of a colloidal aggregation model. The second benchmark considers a network of only bimolecular reactions, where none of the partial propensities are constant. In the third benchmark, we compare PDM and SPDM to SDM on the linear chain model, a weakly coupled reaction network with the minimal degree of coupling, for which SDM was reported to be very efficient \cite{Cao:2004,McCollum:2006}. The fourth benchmark considers the heat-shock response model, a small multiscale (stiff) biological reaction network of fixed size. The benchmark problems are defined in detail in Appendix \ref{app:testproblems}, where also the respective partial-propensity structures $\mathbf{\Pi}$ are given. \\

All tested SSA formulations are implemented in C++ using the random-number generator of the GSL library and compiled using the GNU C++ compiler version 4.0.1 with the O3 optimization flag. All timings are determined using a nanosecond-resolution timer (the {\tt{mach\_absolute\_time()}} system call) on a MacOS X 10.4.11 workstation with a 3\,GHz dual-core Intel Xeon processor, 8\,GB of memory, and a 4\,MB L2 cache. For each test case, we report both the memory requirement and the average CPU time per reaction (i.e.~per time step), $\Theta$. $\Theta$ is defined as the CPU time (identical to wall-clock time in our case) needed to simulate the system up to final time $T$, divided by the total number of reactions executed during the simulation, and averaged over independent runs. The time $\Theta$ does not include the initialization of the data structures (step 1 in Table \ref{our_algorithm}) as this is done only once and is not part of the time loop. \\

We explain the benchmark results in terms of the computational cost of the individual steps of the algorithms. We distinguish three steps: (a) sampling the index of the next reaction, (b) updating the population, and (c) updating the partial propensities (for PDM and SPDM) or the propensities (for SDM). The computational costs of these steps are quantified separately and the overall timings are then explained as a weighted sum of: 
\begin{itemize}
\item \mbox{$\mathcal{C}_\mu$}: The number of operations required to sample the index of the next reaction (for PDM, this is step 2 in Table~\ref{our_algorithm}). 
\item \mbox{$\mathcal{C}_\mathbf{n}$}: The number of elements of the population $\vec{n}$ that need to be updated after executing a reaction (for PDM, this is step 4 in Table~\ref{our_algorithm}). 
\item \mbox{$\mathcal{C}_\r{P}$}: The number of (partial) propensities that need to be updated after executing a reaction (for PDM, this is step 5.2.2 in Table~\ref{our_algorithm}). 
\end{itemize}   
The expressions for these elementary costs are given in Table \ref{table_Cs} as determined by independently fitting models for the scaling of the algorithms to the measured operation counts, averaged over 100 independent runs of each test problem. In all cases, the models used for the computational cost explain the data with a correlation coefficient of at least 0.98. The benchmark results are then explained by fitting the weights of the cost superposition \mbox{$a \mathcal{C}_\mu + b \mathcal{C}_\mathbf{n} + c \mathcal{C}_\r{P}$} to the measured scaling curves $\Theta(N)$ using the expressions given in Table \ref{table_Cs}. In order to preserve the relative weights of the data points, all fits are done on a linear scale, even though the results are plotted on a logarithmic scale for two of the benchmarks. All these fits also have a correlation coefficient of at least 0.98. Explaining the timing results as a superposition of elementary costs allows determining which part of an algorithm is responsible for a particular speedup or scaling behavior, and what the relative contributions of the three algorithmic steps are to the overall computational cost. \\ 

The memory requirements of the algorithms are reported in Table \ref{table_Ms} for all benchmark cases. These numbers were derived analytically from the size of the individual data structures. 

\subsection{Strongly coupled reaction network: colloidal aggregation model}

We use the colloidal aggregation model in Appendix \ref{app:testproblems}\ref{app:aggregation}~\cite{Meakin:1988,Lin:1989,Lin:1990a,Axford:1996,Turner:2005} as a first example of a strongly coupled reaction network. This reaction network can be used to model, e.g., colloidal aggregation of solvated proteins, nano-beads, or viruses. For $N$ chemical species it consists of \mbox{$M = \lfloor{\frac{N^2}{2}}\rfloor$} reactions and the maximum out-degree of the dependency graph is \mbox{$3N-7$} and hence scales with system size (see Table \ref{table_test_cases}).  \\ 

The colloidal aggregation model is simulated up to time \mbox{$T = 100$} with specific probability rates \mbox{$c_{n,m} = 1$} and \mbox{$\bar{c}_{p,q} = 1$}. At time \mbox{$t = 0$}, \mbox{$n_i = N\delta_{1,i}$}. 
The scaling of $\Theta$ for PDM, SPDM, and SDM with system size is shown in Fig.~\ref{Aggregation_model1}(a), averaged over 100 independent runs. \mbox{$\Theta^{\r{PDM}}$} and \mbox{$\Theta^{\r{SPDM}}$} are \mbox{$O(N^{0.5})$} for small $N$ (less than about 100) and \mbox{$O(N)$} for large $N$. \mbox{$\Theta^{\r{SDM}}$} is \mbox{$O(N^{2})$}. The pre-factor of \mbox{$\Theta^{\r{SPDM}}$} is similar to that of \mbox{$\Theta^{\r{PDM}}$}, since in this network \mbox{$\mathcal{C}_\mu$} is not significantly reduced by the dynamic sorting (Table \ref{table_Cs}). The memory requirements of PDM and SPDM are \mbox{$O(N^2)= O(M)$}, that of SDM is \mbox{$O(N^3)= O(NM)$} (Table \ref{table_Ms}). \\  

In summary, the computational costs of both PDM and SPDM are $O(N)$. This scaling is mediated by all three cost components. The use of partial propensities renders the scaling of the sampling cost \mbox{$\mathcal{C}_\mu$} \mbox{$O(N)$} (see Table \ref{table_Cs}). The cost \mbox{$\mathcal{C}_\r{P}$} for updating the partial propensities is \mbox{$O(N^{0.5})$} (Table \ref{table_Cs}), since the use of partial propensities allows formulating a dependency graph over species, rather than reactions, and unimolecular reactions have constant partial propensities. This leads to a smaller number of updates needed as shown in Fig.~\ref{Aggregation_updates}(a). 

\subsection{Strongly coupled network of bimolecular reactions}

The network in Appendix \ref{app:testproblems}\ref{app:tightly} consists of \mbox{$M=\frac{N}{2}(N-1)$} strongly coupled bimolecular reactions, such that none of the partial propensities are constant. Both the minimum and the maximum out-degrees of the dependency graph in this case are \mbox{$4N-10$}, scaling faster with $N$ than in the previous case (see Table \ref{table_test_cases}). \\   

We simulate this network up to time \mbox{$T = 0.001$} with all specific probability rates \mbox{$c_i = 1$}. At \mbox{$t = 0$}, \mbox{$n_i = 100 (\delta_{N-4,i} + \delta_{N-3,i} + \delta_{N-2,i} + \delta_{N-1,i} + \delta_{N,i})$}. The scaling of $\Theta$ for PDM, SPDM, and SDM with system size is shown in Fig.~\ref{Tightly_coupled_model}(b), averaged over 100 independent runs. \mbox{$\Theta^{\r{PDM}}$} and \mbox{$\Theta^{\r{SPDM}}$} are \mbox{$O(N)$}, whereas \mbox{$\Theta^{\r{SDM}}$} is \mbox{$O(N^2)$}. The pre-factors of PDM and SPDM are comparable. The memory requirements of PDM and SPDM are \mbox{$O(N^2)= O(M)$}, that of SDM is \mbox{$O(N^3)= O(NM)$} (see Table \ref{table_Ms}). \\

In summary, the computational costs of PDM and SPDM are $O(N)$ for this strongly coupled, purely bimolecular network. The scaling is again mediated by all three cost components. Grouping the partial propensities renders the sampling cost \mbox{$\mathcal{C}_\mu$} \mbox{$O(N)$} (see Table \ref{table_Cs}). Because none of the partial propensities are constant, the update costs \mbox{$\mathcal{C}_\r{P}$} of PDM and SPDM are $O(N)$, as in SDM, albeit with a pre-factor that is \mbox{$\approx 2.5$} times smaller than that in SDM. One reason for this smaller pre-factor is the smaller number of updates needed upon reactions firing, as shown in Fig.~\ref{Tightly_coupled_updates}(b). This is due to the fact that partial propensities of bimolecular reactions depend on the population of only one species, which reduces the number of combinations that need to be updated.

\subsection{Weakly coupled reaction network: linear chain model}\label{sec:linchain}

We benchmark PDM and SPDM on a weakly coupled model in order to assess their limitations in cases where other SSA formulations might be more efficient. We choose the linear chain model defined in Appendix \ref{app:testproblems}\ref{app:linearchain} since it is the most weakly coupled reaction network possible and it has been used as a model for isolated signal transduction networks~\cite{Albert:2005}. 
For $M$ reactions, it involves the minimum number of species \mbox{$N=M+1$}, and the maximum out-degree of the dependency graph is constant at the minimum possible value of 2 (see Table \ref{table_test_cases}), since every reaction at most influences the population of its only reactant and of the only reactant of the subsequent reaction. \\

We simulate the linear chain model to a final time of \mbox{$T = 1000$} with all specific probability rates \mbox{$c_i = 1$}. At time \mbox{$t = 0$}, \mbox{$n_i = 10000 \delta_{1,i}$}. Figure~\ref{Linear_model1}(c) presents the scaling of the CPU time with system size for PDM, SPDM, and SDM, averaged over 100 independent runs. \mbox{$\Theta^{\r{PDM}}$} scales linearly with $N$ and \mbox{$\Theta^{\r{SPDM}}$} with \mbox{$N^{0.5}$}. \mbox{$\Theta^{\r{SDM}}$} is \mbox{$O(N)$} with a pre-factor that is more than 4 times larger than that of \mbox{$\Theta^{\r{PDM}}$}. This difference in pre-factor is mainly caused by PDM having smaller \mbox{$\mathcal{C}_\mathbf{n}$} and \mbox{$\mathcal{C}_\r{P}$} (Table \ref{table_Cs}). \mbox{$\mathcal{C}_\mu$}, however, scales worse for PDM than for SDM due to the dynamic sorting in SDM. This is overcome in SPDM, where  \mbox{$\mathcal{C}_\mu$} is \mbox{$O(N^{0.5})$}, as in SDM. The memory requirements of SPDM and PDM are \mbox{$O(N)= O(M)$}, that of SDM is \mbox{$O(N^2)= O(NM)$} (Table \ref{table_Ms}). \\

In summary, the computational costs of PDM and SPDM on the weakly coupled linear chain model are governed by (a) updating the population $\vec{n}$ using a sparse stoichiometry representation and (b) never needing to update the partial propensities of unimolecular reactions. Since the linear chain model contains only unimolecular reactions, none of the partial propensities ever needs to be updated, leading to an update cost of \mbox{$\mathcal{C}_\r{P} = 0$} (see Table \ref{table_Cs}). While we have implemented SDM according to the original publication \cite{McCollum:2006}, we note that if one uses a sparse representation of the stoichiometry matrix also in SDM, point (a) vanishes and \mbox{$\mathcal{C}_\mathbf{n}=2$} also for SDM. A sparse-stoichiometry SDM would thus have the same scaling of the computational cost on the linear chain model as would SPDM, outperforming PDM.      

\subsection{Multi-scale biological network: heat-shock response in \textit{Escherichia coli}}

We assess the performance of PDM and SPDM on a small, fixed-size multiscale reaction network. We choose the heat-shock response model since it has also been used to benchmark previous methods, including ODM~\cite{Cao:2004} and SDM~\cite{McCollum:2006}. The model describes one of the mechanisms used by the bacterium {\textit{E.~coli}} to protect itself against a variety of environmental stresses that are potentially harmful to the structural integrity of its proteins. The heat-shock response (HSR) system reacts to this by rapidly synthesizing heat-shock proteins. The heat-shock sigma factor protein \mbox{$\sigma^{32}$} activates the HSR by inducing the transcription of heat-shock genes. The heat-shock response model is a small multiscale reaction network (the specific probability rates span 8 orders of magnitude) with \mbox{$N=28$} chemical species, \mbox{$M=61$} reactions, and a maximum out-degree of the dependency graph of 11 (see Table \ref{table_test_cases}). For a detailed description of the model, we refer to Kurata et al.~\cite{Kurata:2001}\\          

We simulate the HSR model for \mbox{$T = 500$} seconds. During this time, approximately 46 million reactions are executed. 
For a single run, we measure \mbox{$\Theta^{\r{PDM}}=0.256\,\mu$}s and \mbox{$\Theta^{\r{SDM}}=0.272\,\mu$}s. This corresponds to a simulated 3.68 million reactions per second of CPU time for SDM and 3.89 million reactions per second for PDM. Hence, PDM is about 6\% faster than SDM. This speed-up is mainly due to  a smaller \mbox{$\mathcal{C}_\r{P}$} in PDM (see Fig.~\ref{Upsilon_vs_N_plot}(c) for the distribution of updates over all reactions) since the partial propensities of unimolecular reactions never need to be updated. The speed-up is, however, modest because \mbox{$\mathcal{C}_\mu$} of PDM is \mbox{$\approx 4.6$} times larger than that of SDM (Table \ref{table_Cs}). This is due to the fact that 95\% of all reaction firings are caused by a small subset of only 6 reactions. This multiscale network thus strongly benefits from the dynamic sorting used in SDM. This advantage can be recovered in SPDM, where \mbox{$\mathcal{C}_\mu$} is comparable to that of SDM, and \mbox{$\Theta^{\r{SPDM}}=0.245\,\mu$}s (4.08 million reactions per second). This makes SPDM 11\% faster than SDM on this small network.      

\section{Conclusions and Discussion}\label{sec:conclusion}

The stochastic simulation algorithm (SSA)~\cite{Gillespie:1976,Gillespie:1977,Gillespie:1992} is widely used for computational stochastic reaction kinetics in chemistry, physics, biology, and systems biology. It is included in most existing stochastic simulation software packages and is standard in courses on computational chemical kinetics. Due to this importance, several variants of the original SSA formulation have been published that reduce the computational costs of the sampling and update steps. When simulating weakly coupled reaction networks, where the maximum number of reactions that are influenced by any reaction is constant with system size, the computational cost of the sampling step has been reduced to be $O(\log _2 M)$ \cite{Gibson:2000}, where $M$ is the total number of reactions, and even to $O(1)$ under some conditions for the propensity distribution \cite{Slepoy:2008}. Using dependency graphs, also the update step has been reduced to be $O(1)$ for weakly coupled networks \cite{Cao:2004,McCollum:2006,Slepoy:2008}.  
For strongly coupled reaction networks, where the degree of coupling increases with system size and can be as large as the total number of reactions, all previous exact SSA formulations have a computational cost that is $O(M)$. \\

We have introduced a new quantity called {\em partial propensity} and have used it to construct two novel formulations of the exact SSA: PDM and its sorting variant SPDM. Both are algebraically equivalent to DM and yield the same population trajectories $\vec{n}(t)$ as to those produced by DM. In our formulation of partial propensities, we have limited ourselves to elementary chemical reactions. Since their partial propensities depend on the population of at most one species, both new SSA formulations have a computational cost that scales at most linearly with the number of species rather than the number of reactions, independently of the degree of coupling. This is particularly advantageous in strongly coupled reaction networks, where the number of reactions $M$ grows faster than the number of species $N$ with system size. On networks of fixed size, PDM and SPDM are especially efficient when \mbox{$M\gg N$}.
PDM's computational cost is $O(N)$, which is made possible by appropriately grouping the partial propensities in the sampling step and formulating a dependency graph over species rather than reactions in the update step. Moreover, the partial propensities of unimolecular reactions and source reactions are constant and never need to be updated. This further reduces the size of the dependency graph and the computational cost of the update step. To our knowledge, PDM is the first SSA formulation that has a computational cost that is $O(N)$, irrespective of the degree of coupling of the reaction network. In the case of multiscale networks, the computational cost of SPDM is smaller than that of PDM. \\

We have benchmarked PDM and SPDM on four test cases of various degrees of coupling. The first two benchmarks considered strongly coupled networks, where the degree of coupling scales proportionally to the number of species. The third benchmark considered the most weakly coupled network possible, where several other SSA formulations might be more efficient. Finally, the fourth benchmark considered a small biological multiscale network. 
These benchmarks allowed estimating the scaling of the computational cost with system size and the cost contributions from reaction sampling, population update, and partial-propensity update. The results showed that (a) the overall computational costs of PDM and SPDM are $O(N)$, even for strongly coupled networks, (b) on very weakly coupled networks, SPDM is competitive compared to SDM, (c) on multiscale networks SPDM outperforms PDM, and (d) the memory requirements of PDM and SPDM are $O(M)$ in all cases, and hence not larger than those of any other exact SSA formulations. \\

Currently, PDM and SPDM have a number of limitations. The most important limitation is that the presented formulation of partial propensities is only applicable to elementary chemical reactions. Any higher-order chemical reaction can always be broken down into elementary reactions at the expense of increasing system size. In applications such as 
population ecology or social science, the idea of partial propensities can, however, only be used if the (generalized) reactions are at most binary and one species can be factored out, i.e.~if the propensity for every reaction between species $\text{S}_i$ and $\text{S}_j$ can be written as $a_\mu = c_\mu n_i \tilde{h}(n_j)$. Besides this structural limitation, the computational performance of the particular algorithms presented here can be limited in several situations. One of them is the simulation of very small networks, where the overhead of the data structures involved in PDM and SPDM may not be amortized by the gain in efficiency and a simulation using DM may be more efficient. In multiscale networks, where the propensities span several orders of magnitude, PDM is slower than SPDM. In multiscale networks where a small subset ($\ll N$) of all reactions accounts for almost all of the reaction firings, however, the overhead of the data structures involved in SPDM, including their initialization, may not be amortized by the gain in efficiency.
Finally, PDM and SPDM were designed to have a computational cost that scales linearly with the number of species rather than the number of reactions. In reaction networks in which the number of reactions grows sub-linearly with the number of species, this becomes a disadvantage. In such cases, SSA formulations that scale with the number of reactions are favorable.  \\

The classification of reaction networks according to their ``difficulty'' is still largely an open question. Besides system size, degree of coupling, and multiscaling  (spectrum of time scales), there might also be other network properties that influence the computational cost of the various SSA formulations. Automatized selection of the most efficient SSA formulation for a given network would require both a systematic classification of networks and a prediction of the computational cost of SSA formulations based on network properties. 
This might require a more detailed cost analysis of the algorithms and a set of standard benchmark problems that are designed to cover the entire range of performance-relevant parameters. \\

Taken together, our results suggest that PDM and SPDM can potentially offer significant performance improvements especially in strongly coupled networks, including the simulation of colloidal aggregation~\cite{Meakin:1988,Lin:1989,Lin:1990a,Axford:1996,Turner:2005}, Becker-D\"oring-like nucleation-and-growth reactions~\cite{Wattis:2009}, and scale-free biochemical reaction networks, where certain hubs are strongly coupled \cite{Jeong:2000,Strogatz:2001,Albert:2002,Albert:2005}. Finally, the use of partial propensities is not limited to exact SSA formulations, and we also expect approximate methods to benefit from it. 
The software implementations of PDM and SPDM will be made available as open source on the web page of the authors. 

\section{Acknowledgments}
We thank Dr.~Hong Li and Prof.~Dr.~Linda Petzold, University of California at Santa Barbara, for providing the specifications of the heat-shock response model, and the members of the MOSAIC group (ETH Zurich) for fruitful discussions on the manuscript. We also thank the referee for the detailed comments, which greatly helped improving the manuscript, and Jo Helmuth (MOSAIC Group, ETH Zurich) for proofreading. RR thanks Omar Awile for his assistance in optimizing the implementation of PDM. RR was financed by a grant from the Swiss SystemsX.ch initiative, evaluated by the Swiss National Science Foundation. 

