Exact Methods

For a test problem we consider the auto-regulatory network presented in Stochastic Modelling for Systems Biology. There are five species: Gene, P2Gene, Rna, P, and P2, with initial amounts 10, 0, 1, 0, and 0, respectively. There are eight reactions which have mass-action kinetic laws. The table below shows the reactions and propensity factors.


Reaction              Rate constant
----------------------------------------
Gene + P2 → P2Gene    1
P2Gene → Gene + P2    10
Gene → Gene + Rna     0.01
Rna → Rna + P         10
2 P → P2              1
P2 → 2 P              1
Rna → 0               0.1
P → 0                 0.01

Reactions for the auto-regulatory network.

The first figures below shows a single trajectory. A close-up is shown in the next figure. We can see that the system is fairly noisy.


Auto-regulatory system on the time interval [0..50].


Auto-regulatory system on the time interval [0..5].

In order to present a range of problem sizes, we duplicate the species and reactions. For a test problem with 50 species and 80 reactions we have 10 auto-regulatory groups. The reaction propensity factors in each group are scaled by a unit, uniform random deviate. We study systems ranging from 5 to 50,000 species.

The table below shows the performance for various formulations of the direct method. Using a linear search is efficient for a small number of reactions, but does not scale well to larger problems. In the first row we recompute the sum of the propensities at each time step. (This is the original formulation of the direct method.) In the next row we see that immediately updating the sum significantly improves the performance. The following two rows show the effect of ordering the reactions. In the former we periodically sort the reactions and in the latter we swap reactions when modifying the propensities. Ordering the reactions pays off for the largest problem size, but for the rest the overhead outweighs the benefits.

The 2-D search method has the best overall performance. It is fast for small problems and scales well enough to best the more sophisticated methods. Because the auto-regulatory network is so noisy, ordering the reactions hurts the performance of the method.

The binary search on a complete CDF has good performance for the smallest problem size, but has poor scalability. Ordering the reactions is a significant help, but the method is still very slow for large problems. The binary search on a partial, recursive CDF is fairly slow for the smallest problem, but has good scalability. The method is in the running for the second best overall performance.

Because of its complexity, the composition rejection method has poor performance for small problems. However, it has excellent scalability. It edges out the 2-D search method for the test with 80,000 reactions. Although its complexity is independent of the number of reactions, the execution time rises with problem size largely because of caching effects. As with all of the other methods, larger problems and increased storage requirements lead to cache misses. The composition rejection method is tied with the binary search on a partial CDF for the second best overall performance.


Species                                  5   50    500  5,000  50,000 
Reactions                                8   80    800  8,000  80,000
---------------------------------------------------------------------
Linear Search  Delayed update          101  264  1,859 17,145 168,455
Linear Search  Immediate update        109  163    780  6,572  63,113
Linear Search  Complete sort           107  197    976  7,443  22,862
Linear Search  Bubble sort             110  205  1,001  7,420  25,872
2-D Search     Default                 109  130    218    347   1,262
2-D Search     Complete sort           115  148    247    402   1,566
2-D Search     Bubble sort             124  149    220    328   1,674
Binary Search  Complete CDF            105  219  1,196 10,378 103,209
Binary Search  Complete CDF, sorted    114  202    835  3,825  30,273
Binary Search  Partial, recursive CDF  232  328    433    552   1,314
Rejection      Composition             341  365    437    482   1,189

Auto-Regulatory. Direct method. Average time per reaction in nanoseconds.

In the next table we show the performance of the first reaction method. We consider a simple implementation and two implementations that take innovations from the next reaction method. Because a step in the first reaction method has linear computational complexity in the number of reactions, all of the formulations have poor scalability. The simple formulation is fairly slow for small problem sizes. Even for small problems, there is a heavy price for computing the propensity function and an exponential deviate for each reaction. Using the reaction influence graph to reduce recomputing the propensity functions is a moderate help. Storing absolute times instead of the waiting times greatly improves performance. By storing the absolute times, one avoids computing the propensity functions and an exponential deviate for all of the reactions at each time step. Only the reactions influenced by the fired reaction need to be recomputed. However, this formulation is still not competitive with the direct method.


Species               5     50     500    5,000     50,000 
Reactions             8     80     800    8,000     80,000
----------------------------------------------------------
Simple              201  1,968  19,843  159,133  1,789,500
Reaction influence  194  1,510  13,324  110,828    890,948
Absolute time       133    249   1,211   10,368    102,316

Auto-Regulatory. First reaction method. Average time per reaction in nanoseconds.

In the table below we show the performance for various formulations of the next reaction method. Using a linear search is only efficient for a small number of reactions. Manual loop unrolling improves its performance, but it is still not practical for large problems.

The size adaptive and cost adaptive versions of the partition method have pretty good performance. They are competitive with more sophisticated methods up to the test with 800 reactions, but the square root complexity shows in the larger tests.

The binary heap methods have good performance. On 64-bit processors the pair formulation is typically better than the pointer formulation. (Vice-versa for 32-bit processors.)

Using hashing for the priority queue yields the best overall performance for the next reaction method. It is efficient for small problems and has good scalability.


Species                         5   50    500   5,000   50,000 
Reactions                       8   80    800   8,000   80,000
--------------------------------------------------------------
Linear Search  Simple         124  386  2,990  28,902  287,909
Linear Search  Unrolled       120  228  1,116   9,557   94,156
Partition      Fixed size     139  381    582   1,455    5,175
Partition      Size adaptive  163  193    285     500    1,735
Partition      Cost adaptive  124  196    303     537    1,828
Partition      Propensities   146  191    333     723    2,515
Binary Heap    Pointer        166  199    290     413    1,448
Binary Heap    Pair           154  192    272     374    1,304
Hashing        Chaining       151  187    307     320      964

Auto-Regulatory. Next reaction method. Average time per reaction in nanoseconds.

The table below shows the best performing formulation in each category. Only the methods based on a linear search perform poorly. The rest at least offer reasonable performance. The direct method with a 2-D search and the next reaction method that uses a hash table offer the best overall performance. The former is faster up to the test with 800 reactions; the latter has better performance for the large problems.


Species                                         5   50    500  5,000  50,000
Reactions                                       8   80    800  8,000  80,000
----------------------------------------------------------------------------
Direct          Linear search  Complete sort  107  197    976  7,443  22,862
Direct          2-D search     Default        109  130    218    347   1,262
Direct          Binary search  Partial CDF    232  328    433    552   1,314
Direct          Rejection      Composition    341  365    437    482   1,189
First reaction  Linear search  Absolute time  133  249  1,211 10,368 102,316
Next reaction   Linear search  Unrolled       120  228  1,116  9,557  94,156
Next reaction   Partition      Cost adaptive  124  196    303    537   1,828
Next reaction   Binary heap    Pair           154  192    272    374   1,304
Next reaction   Hashing        Chaining       151  187    307    320     964

Auto-Regulatory. Average time per reaction in nanoseconds.

Of course the performance of the various formulations depends upon the problem. The species populations could be highly variable, or fairly stable. The range of propensities could large or small. However, the performance results for the auto-regulatory network are very typical. Most problems give similar results. The biggest difference is that for some systems ordering the reactions is useful when using the direct method. The auto-regulatory system is too noisy for this to improve performance.