Volume 11 Supplement 3
Selected original research articles from the Third International Workshop on Computational Network Biology: Modeling, Analysis, and Control (CNBMAC 2016): systems biology
Stochastic modeling and simulation of reactiondiffusion system with Hill function dynamics
 Minghan Chen†^{1},
 Fei Li†^{1},
 Shuo Wang^{1} and
 Young Cao^{1}Email author
DOI: 10.1186/s1291801704019
© The Author(s) 2017
Published: 14 March 2017
Abstract
Background
Stochastic simulation of reactiondiffusion systems presents great challenges for spatiotemporal biological modeling and simulation. One widely used framework for stochastic simulation of reactiondiffusion systems is reaction diffusion master equation (RDME). Previous studies have discovered that for the RDME, when discretization size approaches zero, reaction time for bimolecular reactions in high dimensional domains tends to infinity.
Results
In this paper, we demonstrate that in the 1D domain, highly nonlinear reaction dynamics given by Hill function may also have dramatic change when discretization size is smaller than a critical value. Moreover, we discuss methods to avoid this problem: smoothing over space, fixed length smoothing over space and a hybrid method.
Conclusion
Our analysis reveals that the switchlike Hill dynamics reduces to a linear function of discretization size when the discretization size is small enough. The three proposed methods could correctly (under certain precision) simulate Hill function dynamics in the microscopic RDME system.
Keywords
Reaction diffusion master equation (RDME) Hill function Stochastic simulation Hybrid methodBackground
Cell reproduction requires elaborate spatial and temporal coordination of crucial events, such as DNA replication, chromosome segregation, and cytokinesis. In cells, protein species are well organized and regulated throughout their life cycles. Theoretical biologists have been using classic chemical reaction rate laws with deterministic ordinary differential equations (ODEs) and partial differential equations (PDEs) to model molecular concentration dynamics in spatiotemporal biological processes. However, wetlab experiments in single cell resolution demonstrate that biological data present considerable variations from cell to cell. The variations arise from the fact that cells are so small that there exist only one or two copies of genes, tens of mRNA molecules and hundreds or thousands of protein molecules [1–3]. At this scale, the traditional way of modeling molecule “concentration” is not applicable. Noise in molecule populations cannot be neglected, as noise may play a significant role in the overall dynamics inside a cell. Therefore, to accurately model the cell cycle mechanism, discrete and stochastic modeling and simulation should be applied.
A convenient strategy to build a stochastic biochemical model is to break a deterministic model into a list of chemical reactions and simulate them with Gillespie’s stochastic simulation algorithm (SSA) [4, 5]. One of the major difficulties in this conversion strategy lies in the propensity calculation of reactions. Gillespie’s SSA is well defined for mass action rate laws. However, in many biochemical models, in addition to mass action rate laws, other phenomenological reaction rate laws are often used. For example, the MichaelisMenten equation [6] and Hill functions [7] are widely used in biological models to model the fast response to signals in regulatory control systems. Although theoretically these phenomenological rate laws may be generated from elementary reactions with mass action rate laws, in practice the detailed mechanisms behind these phenomenological rate laws are not well known and may not be very important. Stochastic modeling and simulation with these phenomenological rate laws are sometimes inevitable.
In recent years, stochastic modeling and simulation for spatiotemporal biological systems, particularly reactiondiffusion systems, have captured more and more attention. Several algorithms and tools [8–11] to model and simulate reactiondiffusion systems have been proposed. These methods can be categorized into two theoretical frameworks: the spatially and temporally continuous Smoluchowski modeling framework [12] and the compartmentbased modeling framework, formulated as the spatially discretized reactiondiffusion master equation (RDME) [13, 14]. The Smoluchowski framework [12, 15, 16] stores the exact position of each molecule and is mathematically fundamental, whereas the RDME is coarsegrained and better suited for large scale simulations [17]. In RDME, the spatial domain is discretized into small compartments. Within each compartment, molecules are considered “wellstirred”. Under the RDME scheme, diffusion is modeled as continuous time random walk on mesh compartments, while reactions fire only among molecules in the same compartment. Stochastic dynamics of the chemical reactions in each compartment is governed by the chemical master equation (CME) [18, 19]. Yet, the CME is computationally impossible to solve for most practical problems. Stochastic simulation methods were then applied to generate realizations of system trajectories. It has been well established that the discretization compartment size for RDME should be smaller than the mean free path of the reactions for the compartment to be considered as wellstirred [20]. In addition, it has been proved that the RDME of bimolecular reactions in 3D domain becomes incorrect and yields unphysical results when the discretization size approaches microscopic scale [21–23].
In this paper, we focus on the stochastic modeling of reactiondiffusion systems with reaction rate laws given by Hill functions. In the Results section, we present our numerical analysis on a toy model of reactiondiffusion system with Hill function dynamics. We will show that the RDME framework of the Hill function dynamics has serious simulation defects when the discretization size approach microscopic limit: When the discretization size is small enough, the typical switching pattern of Hill dynamics becomes linear to the input signal (and the discretization size). Later, we propose potential solutions for the discretization of the reactiondiffusion systems with Hill function rate laws. Finally, we conclude this paper with a discussion on RDME for general nonlinear functions and the hybrid method.
Caulobacter modeling
Caulobacter crescentus captures great interest in the study of asymmetric cell division. When a Caulobacter cell divides, it produces two functionally and morphologically distinct daughter cells. The asymmetric cell division of Caulobacter crescentus requires elaborate temporal and spatial regulations [24–27]. In literature [28–30], four essential “master regulators” of the Caulobacter cell cycle, DnaA, GcrA, CtrA and CcrM, have been identified. These master transcription regulators determine the dynamics of around 200 genes. They oscillate temporally to drive the dynamics of cell cycle. Among them, the molecular mechanisms governing CtrA functions have been well studied. The simulation we are concerned with in this paper is also related to this CtrA module. So we give a brief introduction to it.
In swarmer cells, a twocomponent phosphorelay system (with both CckA and ChpT) phosphorylates the CtrA. Then the chromosomal origin of replication (Cori) is bound by the phosphorylated CtrA (CtrAp) to inhibit the initiation of chromosome replication [31]. Later during the swarmertostalked transition period, CtrAp gets dephosphorylated and degraded, allowing the initiation of chromosome replication again. Thus the CtrA has important impact on the chromosome replication in our model, and should be well regulated.
The regulation of CtrA is achieved by the histidine kinase CckA through the following pathway. An ATPdependent protease, ClpXP, degrades CtrA [32, 33] and is localized to the stalk pole by CpdR. As the nascent stalked cell progresses through the cell cycle, CpdR is phosphorylated by CckA/ChpT, losing it polar localization, and consequently losing its ability to recruit ClpXP protease for CtrA degradation. In addition, CtrA is reactivated through CckA/ChpT phosphorylation [34]. Moreover, the regulatory network of the histidine kinases CckA is influenced by a noncanonical histidine kinase, DivL [35]. DivL promotes CckA kinase, which then phosphorylates and activates CtrA in the swarmer cell. During the swarmertostalked transition period, DivL activity is downregulated, thereby inhibiting CckA kinase activity. As a result, dephosphorylation and degradation of CtrA trigger the initiation of chromosome replication.
In stochastic simulation of the spatiotemporal model of this regulatory network, the phosphorylated CtrA (CtrAp) population switch from a high level in swarmer stage to a low level in stalked stage is not as sharp as expected, shown in Fig. 1 (right). On the other hand, the DivL population level from the stochastic simulation seems similar to that from the deterministic simulation. A simple analysis suggests that the Hill function dynamics, which models the up regulation of CckA kinase activity by DivL, might be the culprit. Further investigation leads to the discovery of the Hill function limitation at small discretization sizes, as analyzed in the next section.
Methods
Reaction diffusion master equation
Before we plunge into Hill functions in reactiondiffusion systems, we will first briefly review mathematical modeling and simulation methods of spatially inhomogeneous stochastic systems.
The dynamics of a spatially inhomogeneous stochastic system has been considered as governed by the reactiondiffusion master equation (RDME), developed in an early work of Gardiner [13]. The RDME framework partitions the spatial domain into small compartments, such that molecules within each compartment can be considered wellstirred. Assume a biochemical system of N species {S _{1},S _{2},…,S _{ N }} and M reactions within a spatial domain Ω, which is partitioned into K grids V _{ k }, k=1,2,…,K. For simplicity, we assume that the space Ω is one dimensional (1D). Each species population, as well as the reactions in the system will have a local copy for each compartment. The state of the reactiondiffusion system at any time t is represented by the vector state vector X(t) = {X _{1,1} (t), X _{1,2}(t), …, X _{1,K }(t),…, X _{ n,k }(t),…,X _{ N,K }(t)}, where X _{ n,k }(t) denotes the molecule population of species S _{ n } in the grid V _{ k } at time t. Reactions in each compartment is governed by the Chemical Master Equation (CME), while diffusion is modeled as random walk across neighboring compartments. Each reaction channel R _{ j } in any compartment k can be characterized by the propensity function a _{ j,k } and the state change vector ν _{ j }≡(ν _{1j },ν _{2j },…,ν _{ Nj }). The dynamics of the diffusion of species S _{ i } from compartment V _{ k } to V _{ j } is formulated by the diffusion propensity function d _{ i,k,j } and the diffusion state change vector μ _{ k,j } similarly. d _{ i,k,j }(x)d t gives the probability that, given X _{ i,k }(t)=x, one molecule of species S _{ i } at grid V _{ k } diffuses into grid V _{ j } in the next infinitesimal time interval [t,t+d t). If j=k±1, then \(d_{i,k,j}(x)=\frac {D}{h^{2}} x\), where D is the diffusion rate coefficient and h is the characteristic length, also called discretization size, of a grid; Otherwise d _{ i,k,j }=0. The state change vector μ _{ k,j } is a vector of length K with −1 in the kth position, 1 in the jth position and 0 everywhere else.
where P(x,tx _{0},t _{0}) denotes the probability that the system state X(t)=x, given that X(t _{0})=x _{0}. The RDME is a set of ODEs that gives one equation for every possible state. It is both theoretically and computationally intractable to solve RDME for practical biochemical systems due to the huge number of possible combinations of states. Instead of solving RDME for the time evolution of the probabilities, we can construct numerical realizations of X(t). A popular method to construct the trajectories of a reactiondiffusion system is to simulate each diffusive jumping and chemical reaction event explicitly. With enough trajectory realizations, we can derive the distribution of each state vector at different times.
The RDME model have been used as an approximation of the Smoluchowski framework in the mesoscopic scale. Furthermore, researches have discovered that in the microscopic limit, bimolecular reactions may be eventually lost when the grid size becomes infinitely small in the three dimensional domain [21, 23]. The RDME framework requires that the two reactant molecules for a bimolecular reaction must be in the same compartment in order to fire a reaction. Intuitively, we may realize that with more discrete compartments, it is less likely for the two molecules to encounter each other at the same compartment in a high dimensional domain. In order to model the reactiondiffusion system with RDME in the microscopic limit, Radek and Chapman [22] derived a formula of meshdependent reaction propensity correction for bimolecular reactions when the discretization size h is larger than a critical size h _{ crit }. This reaction propensity correction formula fails when the discretization size h is smaller than this critical value. Recently, Isaacson [36] proposed a convergent RDME framework (cRDME). In the cRDME framework the diffusion is modeled exactly as in the RDME, while the bimolecular reaction occurs with a nonzero propensity, as long as the distance of the two reactant molecules is less than the reaction radius as defined in the Smoluchowski framework.
In conclusion, the discretization size for the RDME framework should be small enough to avoid discretization error. Yet when the mesh size is less than a critical value, the RDME may become inaccurate for the loss of bimolecular reactions in high dimensional domains. In this paper we will demonstrate that discretization size in space also has great influence on Hill function dynamics in reactiondiffusion systems. The switchlike Hill dynamics breaks even in a 1D domain when the discretization size is small.
Hill function
with V _{ max }=k _{2}[ E]_{0} being the maximum reaction rate and \(K_{m} = \frac {k_{1}+k_{2}}{k_{1}}\) being the Michaelis constant.
with V _{ max } as the maximum reaction rate, K _{ m } as the Michaelis constant, and n as the Hill coefficient. The Hill function is widely used to model “stepregulated” reaction as an activity switch.
Results
The parameters k _{ s }, k _{ d } are the synthesis, degradation rates, respectively, for enzyme species E, and similarly k _{ syn }, k _{ deg } are those for product P. K _{ m } is the Michaelis constant in the Hill function.
where \(\alpha = \frac {k_{s}}{k_{d}}L\) denotes the mean value of the total number of enzyme E molecules in the domain.
The other probability terms are not important in the analysis.
Notice that α/K is the mean population of enzyme E in the ith bin. The Hill function of the product P synthesis is now reduced to a linear function of the enzyme E population in the ith bin.
The loglog plot (Fig. 3, right) shows that when the discretization size is small enough, the total product P population is a linear function of discretization size. The slope of the loglog plot is about 1.0 at small discretization size h, regardless of K _{ m }.
Moreover, simulation results show that when the mean enzyme E population is less than the constant K _{ m } in the Hill function (K _{ m }>α), the population of product P increases slightly before the Hill function dynamics breaks at small discretization sizes. Note that the Hill function dynamics show a concave shape with respect to enzyme E population when the enzyme E population is smaller than the Michaelis constant K _{ m }. Therefore, it is reasonable that the product P population in this reactiondiffusion model increases slightly when the Michaelis constant K _{ m } is larger than the mean enzyme E population α.
the Hill dynamics reduce to a linear function.
Equivalently, in order for the Hill function dynamics to work well, the discretization number K should be less than or equal to this threshold. However, the coarse discretization from a small K leads to spatial error. Two potential solutions to this discretization dilemma are proposed next.
Discussion
From the previous analysis, the Hill dynamics in RDME systems reduces to a linear function due to the lack of intermediate states — the discrete population in each individual bin yields an integer value (0 or 1) for the Hill function. Thus a natural solution to it is to generate intermediate states by a smoothing technique that averages the population over neighboring bins when calculating the reaction propensity.
To model a RDME system in high dimensions with fine discretization, previous studies [21] have suggested relaxing the samecompartment reaction assumption and allowing reactions within neighboring compartments. The next subsection shows that allowing reactions within neighboring compartments is equivalent to smoothing over neighboring compartments.
Smooth over neighboring bins
A natural technique that bridges the discrete and continuous models is to smooth the spatial population by taking the average of neighboring bins. Consider first smoothing the enzyme E population within the neighboring m bins (including the bin itself) when calculating the reaction propensity.
As before, only the term \({P}_{E}^{(i)}(0; m)\) is important.
Convergent hill function dynamics in reactiondiffusion systems
The previous subsection demonstrates that a sufficiently small discretization size h will still break the Hill dynamics even with the strategy of smoothing over a fixed number of bins, thus the number of bins needs to vary with the discretization size.
Inspired by the convergentRDME framework [36], a remedy for the failure of Hill function dynamics in reactiondiffusion systems is to smooth the population over bins within a certain distance.
Following the terminology in the convergentRDME framework [36], the “reaction radius ρ” of the Hill function dynamics is defined as ρ=m·h, where m is given in (23).
Conclusions
Motivated by the misbehavior of DivLCckA dynamics in the stochastic simulation of the Caulobacter crescentus cell cycle model, a study of the Hill function dynamics in reactiondiffusion systems reveals that when the discretization size is small enough, the switchlike behavior of Hill function dynamics reduces to a linear function of input signal and discretization size. A proposed fixed length smoothing method, which allows chemical reactions to occur with reactant molecules within a distance of fixed length, the “reaction radius”of the Hill function dynamics, seems to give a very good remedy to this problem.
It is known that in high dimensions bimolecular reactions are lost with the RDME in the microscopic limit [21]. This work shows that onedimensional Hill function dynamics in a RDME framework gives a similar challenge when the discretization size is small enough. The conjecture is that the problem lies in the RDME requirement that reactions only fire with the reactant molecules in the same discrete compartment.
Furthermore, this defect in RDME at the microscopic limit is believed to be a common scenario for all highly nonlinear reaction dynamics. Theoretical biologists have developed many highly nonlinear reaction dynamics that need special attention when converted to stochastic models.
where P(·) is the probability under the distribution when the virtual fast system (25) converges to stochastic partial equilibrium [37].
(28) highlights the mismatch between the RDME framework and the deterministic model.
Hybrid method

Two uniform random numbers, r _{1} and r _{2} in U(0,1), are generated.

Solve the ODE system for S _{ fast } and find the root τ for the integral equation:$$ \int^{t+\tau}_{t} a_{tot}(\mathbf{x},s)ds+\log(r_{1}) = 0, $$(29)
where a _{ tot }(x,t) is the sum of propensities of all reactions in S _{ slow }. Because x varies with t in the ODE system, a _{ tot }(x,t) is a function of t as well.

μ is selected as the smallest integer satisfying$$ \sum^{\mu}_{i=1}a_{i}(\mathbf{x},t) > r_{2} a_{tot}(\mathbf{x},t). $$(30)

Update \(\mathbf {x}\leftarrow \mathbf {x}+\mathbf {\nu }_{\mu }\).

Return to step 1) if stopping condition is not reached.
where we note that log(r _{1}) is negative and a _{ tot } is always nonnegative. In the hybrid simulation, for each step we start from the current time t and numerically [39] integrate the original ODEs (31) and the extra integral Eq. (32). The integration stops when z(t+τ)=0. As a result, τ is the solution to (29). This procedure can be numerically simulated using standard ODE solvers combined with rootfinding functions, such as the LSODAR [39]. Note that since z is an integration variable, one may choose to omit it from the error control mechanism [40]. Adding this extra variable will not greatly affect the efficiency.
Numerical results certainly suggest that the hybrid method has great potential in stochastic simulation of RD systems. We would like to note that great details still need to be studied, but that is not the focus for this paper.
Abbreviations
 CME:

Chemical master equation
 ODE:

Ordinary differential equation
 PDE:

Partial differential equation
 RDME:

Reaction diffusion master equation
 SSA:

Stochastic simulation algorithm
Declarations
Funding
This work was partially supported by the National Science Foundation award DMS1225160, CCF0953590, CCF1526666, and MCB1613741. In particular, the publication of this paper is directly funded from CCF1526666, and MCB1613741.
Availability of data and materials
Not applicable.
Authors’ contributions
FL initially realized the simulation error described in this paper. FL, MC and YC then designed the toy model and analyzed the numerical error caused by Hill function simulation. Later SW joined to help the implementation of hybrid simulation. FL and MC together drafted the manuscript, and YC gave critical revisions on the writing. All authors have read and approved the final manuscript.
Competing interests
The authors declare that they have no competing interests.
Consent for publication
All authors consent to publish this work through BMC’17.
Ethics approval and consent to participate
Not applicable.
About this supplement
This article has been published as part of BMC Systems Biology Volume 11 Supplement 3, 2017: Selected original research articles from the Third International Workshop on Computational Network Biology: Modeling, Analysis, and Control (CNBMAC 2016): systems biology. The full contents of the supplement are available online at http://bmcsystbiol.biomedcentral.com/articles/supplements/volume11supplement3.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License(http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver(http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
Authors’ Affiliations
References
 McAdams H, Arkin A. Stochastic mechanisms in gene expression. Proc Natl Acad Sci. 1997; 94(3):814–9. http://www.pnas.org/content/94/3/814.
 Fedoroff N, Fontana W. Small numbers of big molecules. Science. 2002; 297(5584):1129–31. http://www.sciencemag.org/content/297/5584/1129.
 Samoilov M, Plyasunov S, Arkin AP. Stochastic amplification and signaling in enzymatic futile cycles through noiseinduced bistability with oscillations. Proc Natl Acad Sci U S A. 2005; 102(7):2310–5. http://www.pnas.org/content/102/7/2310.
 Gillespie DT. A general method for numerically simulating the stochastic time evolution of coupled chemical reactions. J Comput Phys. 1976; 22(4):403–34.View ArticleGoogle Scholar
 Gillespie DT. Exact stochastic simulation of coupled chemical reactions. J Phys Chem. 1977; 81(25):2340–61.View ArticleGoogle Scholar
 Michaelis L, Menten ML. Die Kinetik der Invertinwirkung. Biochem Z. 1913; 49:333–69.Google Scholar
 Hill AV. The possible effects of the aggregation of the molecules of haemoglobin on its dissociation curves. J Physiol. 1910; 40(Suppl):iv–vii.Google Scholar
 Hattne J, Fange D, Elf J. Stochastic reactiondiffusion simulation with MesoRD. Bioinformatics. 2005; 21(12):2923–4.View ArticlePubMedGoogle Scholar
 Andrews SS, Bray D. Stochastic simulation of chemical reactions with spatial resolution and single molecule detail. Phys Biol. 2004; 1:137–51.View ArticlePubMedGoogle Scholar
 van Zon JS, ten Wolde PR. Green’sfunction reaction dynamics: A particlebased approach for simulating biochemical networks in time and space. J Chem Phys. 2005; 123(23):234910. http://scitation.aip.org/content/aip/journal/jcp/123/23/10.1063/1.2137716.
 Novère NL, Shimizu TS. StochSim: modelling of stochastic biomolecular processes. Bioinformatics. 2001; 17:575–6.View ArticlePubMedGoogle Scholar
 von Smoluchowski M. Zur kinetischen Theorie der Brownschen Molekularbewegung und der Suspensionen. Annalen der Physik. 1906; 326(14):756–80. http://dx.doi.org/10.1002/andp.19063261405.
 Gardiner CW, McNeil KJ, Walls DF, Matheson IS. Correlations in stochastic theories of chemical reactions. J Stat Phys. 1976; 14:307–31.View ArticleGoogle Scholar
 Nicolis G, Prigogine I. Selforganization in nonequilibrium systems : from dissipative structures to order through fluctuations. New York: A WileyInterscience Publication; 1977. http://opac.inria.fr/record=b1078628.Google Scholar
 doi M. Stochastic theory of diffusioncontrolled reaction. J Phys A Math General. 1976; 9(9):1479. http://stacks.iop.org/03054470/9/i=9/a=009.
 Keizer J. Nonequilibrium statistical thermodynamics and the effect of diffusion on chemical reaction rates. J Phys Chem. 1982; 86(26):5052–67. http://dx.doi.org/10.1021/j100223a004.
 Fange D, Berg OG, Sjöberg P, Elf J. Stochastic reactiondiffusion kinetics in the microscopic limit. Proc Natl Acad Sci. 2010; 107(46):19820–5. http://www.pnas.org/content/107/46/19820.
 McQuarrie DA. Stochastic approach to chemical kinetics. J Appl Probab. 1967; 4(3):413–78. http://www.jstor.org/stable/3212214.
 Gillespie DT. A rigorous derivation of the chemical master equation. Physica A: Stat Mech Appl. 1992; 188(1–3):404–25.View ArticleGoogle Scholar
 Baras F, Mansour MM. Reactiondiffusion master equation: a comparison with microscopic simulations. Phys Rev E. 1996; 54:6139–48. http://link.aps.org/doi/10.1103/PhysRevE.54.6139.
 Isaacson SA. The reactiondiffusion master equation as an asymptotic approximation of diffusion to a small target. SIAM J Appl Math. 2009; 70(1):77–111.View ArticleGoogle Scholar
 Erban R, Chapman SJ. Stochastic modelling of reactiondiffusion processes: algorithms for bimolecular reactions. Phys Biol. 2009; 6(4):046001. http://stacks.iop.org/14783975/6/i=4/a=046001.
 Hellander S, Hellander A, Petzold L. Reactiondiffusion master equation in the microscopic limit. Phys Rev E. 2012; 85:042901. http://link.aps.org/doi/10.1103/PhysRevE.85.042901.
 Li S, Brazhnik P, Sobral B, Tyson JJ. A Quantitative Study of the Division Cycle of Caulobacter crescentus, Stalked Cells. PLoS Comput Biol. 2008; 4(1):e9. http://dx.plos.org/10.1371%252Fjournal.pcbi.0040009.
 Li S, Brazhnik P, Sobral B, Tyson JJ. Temporal controls of the asymmetric cell division cycle in Caulobacter crescentus. PLoS Comput Biol. 2009; 5(8):000463. http://dx.doi.org/10.1371%252Fjournal.pcbi.1000463.
 Subramanian K, Paul MR, Tyson JJ. Potential role of a bistable histidine kinase switch in the asymmetric division cycle of Caulobacter crescentus. PLoS Comput Biol. 2013; 9(9):003221. http://dx.doi.org/10.1371%2Fjournal.pcbi.1003221.
 Subramanian K, Paul MR, Tyson JJ. Dynamical localization of DivL and PleC in the asymmetric Division cycle of Caulobacter crescentus: a theoretical investigation of alternative models. PLoS Comput Biol. 2015; 11(7):004348.View ArticleGoogle Scholar
 Collier J, Murray SR, Shapiro L. DnaA couples DNA replication and the expression of two cell cycle master regulators. EMBO J. 2006; 25(2):346–56.View ArticlePubMedPubMed CentralGoogle Scholar
 Collier J, Shapiro L. Spatial complexity and control of a bacterial cell cycle. Curr Opin Biotechnol. 2007; 18(4):333–40. http://www.sciencedirect.com/science/article/pii/S0958166907000894.View ArticlePubMedPubMed CentralGoogle Scholar
 Holtzendorff J, Hung D, Brende P, Reisenauer A, Viollier PH, McAdams HH, et al. Oscillating global regulators control the genetic circuit driving a bacterial cell cycle. Science. 2004; 304(5673):983–7. http://www.sciencemag.org/content/304/5673/983.View ArticlePubMedGoogle Scholar
 Quon KC, Yang B, Domian IJ, Shapiro L, Marczynski GT. Negative control of bacterial DNA replication by a cell cycle regulatory protein that binds at the chromosome origin. Proc Natl Acad Sci. 1998; 95(1):120–5. http://www.pnas.org/content/95/1/120.View ArticlePubMedPubMed CentralGoogle Scholar
 McGrath PT, Iniesta AA, Ryan KR, Shapiro L, McAdams HH. A dynamically localized protease complex and a polar specificity factor control a cell cycle master regulator. Cell. 2006; 124(3):535–47. http://www.sciencedirect.com/science/article/pii/S0092867406000663.View ArticlePubMedGoogle Scholar
 Jenal U, Fuchs T. An essential protease involved in bacterial cellcycle control. EMBO J. 1998; 17(19):5658–69.View ArticlePubMedPubMed CentralGoogle Scholar
 Iniesta AA, McGrath PT, Reisenauer A, McAdams HH, Shapiro L. A phosphosignaling pathway controls the localization and activity of a protease complex critical for bacterial cell cycle progression. Proc Natl Acad Sci. 2006; 103(29):10935–40. http://www.pnas.org/content/103/29/10935.View ArticlePubMedPubMed CentralGoogle Scholar
 Tsokos C, Perchuk B, Laub M. A dynamic complex of signaling proteins uses polar localization to regulate cellfate asymmetry in Caulobacter crescentus. Developmental Cell. 2011; 20(3):329–41.View ArticlePubMedPubMed CentralGoogle Scholar
 Isaacson SA. A convergent reactiondiffusion master equation. J Chem Phys. 2013; 139(5):054101. http://scitation.aip.org/content/aip/journal/jcp/139/5/10.1063/1.4816377.View ArticlePubMedGoogle Scholar
 Cao Y, Gillespie DT, Petzold LR. The slowscale stochastic simulation algorithm. J Chem Phys. 2005; 122(1):014116.View ArticleGoogle Scholar
 Haseltine EL, Rawlings JB. Approximate simulation of coupled fast and slow reactions for stochastic chemical kinetics. J Chem Phys. 2002; 117(15):6959–69.View ArticleGoogle Scholar
 Hindmarsh AC. ODEPACK, a systematized collection of ODE solvers. IMACS Trans Sci Comput Amsterdam. 1983; 1:55–64.Google Scholar
 Petzold LR. A Description of DASSL: a differential/algebraic system solver, proceeding of the 1st IMACS World Congress. Montreal. 1982;1:6568.