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.