Modeling suggests that gene circuit architecture controls phenotypic variability in a bacterial persistence network
© Koh and Dunlop; licensee BioMed Central Ltd. 2012
Received: 19 December 2011
Accepted: 20 May 2012
Published: 20 May 2012
Bacterial persistence is a non-inherited bet-hedging mechanism where a subpopulation of cells enters a dormant state, allowing those cells to survive environmental stress such as treatment with antibiotics. Persister cells are not mutants; they are formed by natural stochastic variation in gene expression. Understanding how regulatory architecture influences the level of phenotypic variability can help us explain how the frequency of persistence events can be tuned.
We present a model of the regulatory network controlling the HipBA toxin-antitoxin system from Escherichia coli. Using a biologically realistic model we first determine that the persistence phenotype is not the result of bistability within the network. Next, we develop a stochastic model and show that cells can enter persistence due to random fluctuations in transcription, translation, degradation, and complex formation. We then examine alternative gene circuit architectures for controlling hipBA expression and show that networks with more noise (more persisters) and less noise (fewer persisters) are straightforward to achieve. Thus, we propose that the gene circuit architecture can be used to tune the frequency of persistence, a trait that can be selected for by evolution.
We develop deterministic and stochastic models describing how the regulation of toxin and antitoxin expression influences phenotypic variation within a population. Persistence events are the result of stochastic fluctuations in toxin levels that cross a threshold, and their frequency is controlled by the regulatory topology governing gene expression.
KeywordsPersister Toxin-antitoxin Gene regulatory network Feedback
Gene expression is controlled by regulatory networks that influence the mean levels, dynamics, and noise distributions of proteins expressed within a single cell. The outputs of these networks are under selective pressure; thus a regulatory architecture that results in beneficial traits can be selected for by evolution. A key question in systems biology is how the architecture of a gene regulatory network influences the dynamics of gene expression. This question has been explored extensively using mathematical modeling . However, a subtler question is how the architecture of a gene circuit influences the variability in gene expression, and what the implications are for population fitness. Previous studies have shown that similar gene circuit architectures can produce vastly different noise profiles [2, 3]. It is clear from systems-level studies that noise in gene networks can be controlled, or selected for, by evolution [4–7]. For example, stress-response genes in Saccharomyces cerevisiae have been shown to exhibit significant stochastic variation ; similar results have been found for genes that respond to environmental changes . Furthermore, genes that are lethal when deleted exhibit much lower than average levels of stochastic variation . The regulatory architecture of a network also plays an important role in controlling noise [1–4, 11–15]. This has been shown specifically for drug-induced stress through several studies demonstrating that increased phenotypic variability can provide a selective advantage [16–19]. The regulation of noise can have dramatic implications when it controls physiological processes important for stress response and survival.
Here, we explore the role of network architecture on the noise properties of a regulatory circuit controlling bacterial persistence. Persistence is a non-inherited mechanism by which bacteria tolerate environmental stress, such as treatment by antibiotics. Cells are able to stochastically switch between a dormant state known as persistence and a regular growth state. In the persistence state, cell growth slows drastically and the cell is therefore immune to treatment by antimicrobial agents that target growth. Examples include beta-lactams, which interfere with cell wall biosynthesis and aminoglycosides, which interrupt translation . Because the cells are in a dormant state, they can effectively evade the drugs. When the cells switch back to the regular growth state, they resume replication and are no longer tolerant to antimicrobial agents. It is important to understand that persistence is not a genetic change; both persisters and normal cells have identical genetic code. Instead, persistence is a transient state that a small subset of the population enters due to phenotypic variation. By maintaining subpopulations of normal and persister cells, the whole population hedges against unlikely but catastrophic events, while still maintaining near optimal growth at the population level.
Persistence plays an important role in chronic infections. High-persistence (or hip) mutants are found in Pseudomonas aeruginosa, Candida albicans, Escherichia coli, and Mycobacterium tuberculosis[21, 22]. At any given time, a typical bacterial population will have between 10-7 and 10-5 cells in the persistence state . High persistence mutants result in 10-4 to 10-2 portion of cells in the persistence state . This suggests that although it is possible to increase the number of persisters through mutations, their level has been tuned to balance tolerance to environmental threats with a maximal growth rate.
Toxin-antitoxin modules play a key role in the formation of persister cells. Previous work has shown that toxins can induce the dormant state by inhibiting important cellular processes, most commonly mRNA translation . Though toxins have historically been associated with programmed cell death, recent studies show that toxin-antitoxin loci can function to moderate global levels of translation and replication in cells that must survive in stressful environments . For example, the toxins relE and mazF cleave mRNA, thus inhibiting translation and cell replication. This helps the cell cope with nutritional stress; if there is a temporary shortage of nutrients, some of the population will survive. The action of the toxin can be countered by the presence of a small antitoxin molecule which binds to the large toxin protein, disabling its function . The toxin is stable relative to the antitoxin, whose rapid degradation and production ensure a fast dynamic response. This creates a two-state system, where an excess of toxin will induce dormancy and persistence, but enough antitoxin keeps the cell in its normal growth state. Many pathogenic bacteria have multiple toxin-antitoxin loci. For example, M. tuberculosis has 88 known toxin-antitoxin systems . There may be interactions between these different toxin-antitoxin loci that further complicate our current understanding of persistence.
The regulatory architecture of toxin-antitoxin systems is highly conserved across bacterial species . In most cases, the antitoxin gene precedes the toxin gene, with the two loci expressed on the same operon. Because of this, the two genes share a regulatory structure, which can serve to correlate noise in toxin and antitoxin expression. In addition, the antitoxin protein is often a transcription factor, which autorepresses the promoter that controls expression of the two genes. This negative feedback-based inhibition results in low levels of protein expression, which increases the potential for noise in the system. Without feedback, higher toxin and antitoxin expression would mitigate the effect of noise. This tradeoff in noise, which is reduced due to coupled transcription, but elevated due to negative feedback, is the focus of the present study.
The HipBA toxin-antitoxin system from E. coli is a well-known persistence mechanism. The HipA toxin causes cell stasis, but HipB can inactivate HipA and create a non-toxic complex. A HipB dimer binds to two copies of HipA, rendering it neutral through conformational inactivation and sequestration . Here, we use the HipBA toxin-antitoxin system as a model for understanding the dynamics and regulatory processes that govern bacterial persistence. HipA is highly conserved among Gram-negative bacteria , suggesting that phenotypic variation in toxin-antitoxin expression is a common mechanism for persistence.
The aim of this study is to identify the specific gene expression dynamics that govern persistence. In particular, we ask how the regulatory architecture of the gene circuit leads to noise, and whether this noise is subject to evolutionary tuning. Knowledge of how, why, and when cells switch to persistence can help guide studies on treatment strategies to reduce or eliminate the number of cells that enter persistence.
Chemical reaction parameters for the HipBA system
M + P
M + P'
M + P''
HipB binds to promoter
P + B2
approximated based on 
P + B2
approximated based on 
HipB-HipA complex binds to promoter
P + AB2A
P + AB2A
M + B
M + A
B + B
assumed monomeric form to be uncommon
HipB dimer degradation
HipB-HipA complex association
B2 + 2A
approximated based on 
HipB-HipA complex dissociation
B2 + 2A
approximated based on 
HipB-HipA complex degradation
A represents the HipA protein; B is HipB. B2 is the dimerized form of HipB. AB2A is the HipA-HipB toxin-antitoxin complex. M is the mRNA transcript from hipBA. P is the promoter of hipBA with no proteins bound, P’ has B2 bound, and P” has AB2A bound. The equations for the rate of change of P, P’, and P” describe how the promoter switches between states with nothing, B2, and AB2A bound. mRNA is transcribed from all three promoter states at different rates and is also degraded. HipA is translated from mRNA and degraded by a protease. HipB is also translated from mRNA and subsequently binds to a second copy of itself to form the HipB dimer, which can bind to and repress the promoter. Dimerization between two HipB molecules is modeled as irreversible because unbinding is slow relative to the binding rate [27, 29]. Two copies of HipA bind to one copy of the HipB dimer to form the HipB-HipA complex, which can bind to and represses the promoter.
To analyze the possibility of multiple steady state solutions, we developed a reduced order model and conducted phase plane analysis. The dynamics of the promoter states, mRNA, B, and B2 are fast relative to A and AB2A. Thus, we assumed that the other states were at equilibrium, and therefore time derivatives equal to zero. The steady state concentrations of the fast variables were used in a two-dimensional model that describes the rate of change of A and AB2A. We next used a phase portrait to plot the nullclines, which are the lines where d[A]/dt = 0 and d[AB2A]/dt = 0. The points where the nullclines cross are the equilibrium points of the system. If the lines cross more than once, multiple equilibrium solutions are possible. The full equations for the reduced order model are given in the Additional file 1.
Next, we conducted two parametric studies to verify that the system dynamics are monostable for a broad range of biologically realistic parameter values. First, we varied single parameters and checked for the existence of multiple stable states. For each of the parameters in the model we chose that parameter from a log normal distribution using the range given in Table 1. All other parameters were held constant. Using this set of parameters, the system was simulated for 500 distinct initial conditions, which were generated randomly with initial mRNA concentration varying uniformly between 0 and 100 and all proteins varying uniformly between 0 and 1000. Promoter concentrations were held constant for all simulations. We then checked, numerically, whether starting the system at different initial conditions generated any solutions that were more than 1% away from any others. This test was repeated 100 times for each of the 17 model parameters. We found no solutions that varied more than 1% for different initial conditions, suggesting that only monostable solutions exist.
Next, we allowed all system parameters to vary simultaneously. Specifically, all parameters were selected from log normal distributions using the ranges given in Table 1. Initial conditions were selected as described in the single parameter study. We tested 1000 different combinations of parameters and checked for the existence of multiple stable states using the 1% metric described above. For all parameter combinations the systems converged to a single equilibrium point.
Cells are defined as entering persistence when the number of free HipA toxins exceeds the number of free HipB antitoxins. We used the ratio of free HipA molecules in the cell divided by the sum of free HipA and free HipB molecules to quantify entry into persistence. Hereafter, we will refer to this quantity as R. When R exceeds 0.5 the cell is a persister; below this value the cell is in the normal growth state. This threshold-based approach is consistent with experimental findings from , where the levels of HipA and HipB were controlled independently using a plasmid expression system and were found to depend on the concentration of HipB.
The “uncoupled transcription” model replaces the original three promoter states P, P’, and P” with six promoter states, three for hipB and three for hipA. In addition, the mRNA transcripts MA and MB are now two separate states in the model. Thus, the new model has 12 state variables, but uses the same reaction constants as the native system for the purpose of a controlled comparison (Additional file 1).
The “no feedback” model removes repression of the hipBA promoter by B2 and AB2A. To model this we considered only the P promoter state, eliminating the P’ and P” states from the model. All other reactions and reaction rates are the same as in the original system (Additional file 1).
Results and discussion
In order to analyze the dynamics of persister formation in the HipBA system we developed a mathematical model based on the regulatory architecture known to control HipB and HipA expression. We first asked how the dynamics of the system led to the formation of persister cells. Next, we studied alternative network architectures to quantify how entry into persistence depends upon gene regulatory structure. In order to address these questions, we developed a biologically realistic model. The system explicitly models promoter states, the binding and unbinding of transcription factors, transcription, translation, complex formation, and degradation. In contrast to previous models [25, 31], we consider dimerization, complex formation, multiple modes of repression, and active degradation of the toxin and antitoxins; these processes are modeled based on the physiological findings from experimental studies (Methods).
We first asked how the HipBA regulatory architecture achieves distinct subpopulations of persister and normal cells. A potential mechanism for generating two populations within a group of cells is bistability. There is experimental evidence that isogenic populations can generate bimodal distributions to allow for phenotypic diversity [34, 35]. This strategy is beneficial when only a subset of the cells needs to express a particular mechanism, but those cells need to be fully committed to their fate. Positive feedback is known to generate bistable states and can arise from a double negative feedback loop . Therefore, a potential function of the HipBA regulatory network could be to generate two stable states through the use of two negative feedback loops, which act in combination as a positive feedback loop (Figure 1b). In principle, repression of the promoter by the HipB dimer or HipB-HipA complex could lead to a build up of the HipA toxin, and consequently persistence, because the half-life of HipA exceeds that of HipB [21, 36]. Alternatively, the higher translation rate of HipB could lead to an excess of antitoxins, leading to the normal growth state. Stochastic fluctuations in gene expression could cause the system to switch between these two states. A previous study proposed a model for persistence based on high cooperativity in a Hill function as a mechanism that generates bistable dynamics .
Using our detailed mechanistic model of the biochemical reactions governing HipB and HipA expression we found that the system was monostable for biologically realistic parameter ranges, therefore bistability is not the source of co-existing persister and normal cells. To check for bistable dynamics, we first used time scale separation to develop a reduced order model (Methods, Additional file 1). The dynamics of HipA (A) and the HipB-HipA complex (AB2A) were slow relative to the other states in the system. Thus, we developed a reduced order model that assumed other chemical reactants were at steady state relative to A and AB2A. We then plotted the nullclines for A and AB2A on a phase portrait and showed that, for realistic parameter ranges, they intersect only once (Figure 1c). This single intersection point indicates that only one equilibrium solution exists, thus the system is not bistable.
In order to rule out the possibility that the absence of bistability was the result of the specific parameters used in the model, we conducted two parametric studies (Methods). First, we varied single parameters within a biologically realistic range (Table 1) and tested for bistable behavior over a broad range of initial conditions. In all cases, the solutions converged to a single monostable equilibrium point. Next, we allowed all system parameters to vary at once, and simulated many possible combinations of parameters. Again, solutions for all parameters converged to a single stable point. Through a combination of reduced order system analysis and parametric studies, we find no evidence of bistability in our model of the HipBA system.
An alternative mechanism by which cells can enter persistence is through stochastic fluctuations in gene expression. Random noise in the expression of HipB and HipA can generate phenotypic variability within the population. By chance, some cells within the population will have an excess of the toxin relative to the antitoxin and will enter persistence. To explore the role of phenotypic variability in persister formation, we developed a stochastic model based on the chemical reactions used in the deterministic model. The probabilistic nature of this model more accurately represents the natural fluctuations in the HipBA system.
In order for the HipA toxin to be effective, an individual cell would have to have an excess of free HipA toxins relative to the number of free HipB antitoxins. Thus, the ratio of free HipA molecules to the total number of free HipA and HipB molecules, which we define as R, sets a threshold for persistence. When R exceeds 0.5 a cell has an excess of toxin and can enter persistence. Recent experimental findings suggest that a threshold-based mechanism for persistence, as opposed to bistability, is an accurate representation of the biological origins of persistence . The authors showed that the time spent in the persistence state was proportional to the concentration of excess HipA. Our study examines the entry into persistence, however the duration of the growth arrest period is not calculated by the model, as the dynamics are only valid for non-persister cells. Therefore, our model can be used to simulate the distributions of HipA and HipB and this information can be used to calculate entry into persistence, but not the duration of the growth arrest state.
A system with increased persistence would be better suited for conditions where extreme environmental stress occurs frequently or for extended periods of time. Although the population growth rate would be severely compromised, cells would have an increased likelihood of surviving extreme or long-term environmental stresses, such as long-term nutrient deprivation or antibiotic treatment. Conversely, a system with decreased persistence would benefit from increased growth rates and thrive in environments where stresses are few and far between. A previous model of persistence has shown that the optimal frequency of persistence events is closely tied to the frequency of environmental change . The regulatory topology of the HipBA network sets a frequency of entry into persistence, which may be a strong indicator of the frequency with which adverse environments are encountered.
Evolutionarily, it would be possible to achieve either of the alternative circuit topologies discussed here through straightforward mutation or duplication events. Given that stochastic fluctuations in phenotypic states are the likely source of persisters, it is necessary for the regulatory architecture to produce sufficient variability to insure against rare but catastrophic environmental stresses. This suggests that the HipBA toxin-antitoxin system has evolved to allow a specific amount of noise, and thus persistence, to balance between optimal growth and survival against environmental threats.
We have developed a model of the regulatory interactions that control expression of the hipBA operon. The model incorporates recent experimental findings on toxin-antitoxin complex formation , includes active degradation of proteins [21, 31, 32], and repression of hipBA gene expression by both the HipB dimer and the toxin-antitoxin complex . Using a deterministic model based on the biochemical reactions we find that the system exhibits monostable behavior for biologically realistic parameters. Thus, stochastic fluctuations in expression of HipB and HipA are likely responsible for the spontaneous formation of persisters. A stochastic model of the biochemical system demonstrates that persister formation is possible when R, the ratio of free HipA to free HipA plus free HipB exceeds a threshold. By comparing two alternative gene circuit architectures we demonstrate that the level of noise, and thus the frequency of persistence, is highly dependent on the regulatory topology. Through mutations it would be possible to achieve systems with either higher or lower noise than the wild type system. Populations must balance the tradeoff between frequency of persistence and growth. We conclude that straightforward changes to the regulatory architecture and associated parameters influence noise levels, which define the frequency of persistence, suggesting that the wild type persistence frequency can be the subject of evolutionary tuning.
We thank all group members and members of the Hill group for useful discussions. This work was supported by funding from the University of Vermont to MJD.
- Alon U: An Introduction to Systems Biology: Design Principles of Biological Circuits.2007.Google Scholar
- Cağatay T, et al., et al.: Architecture-dependent noise discriminates functionally analogous differentiation circuits. Cell 2009,139(3):512-522. 10.1016/j.cell.2009.07.046View ArticleGoogle Scholar
- Kittisopikul M, Suel GM: Biological role of noise encoded in a genetic network motif. Proc Natl Acad Sci U S A 2010,107(30):13300-13305. 10.1073/pnas.1003975107View ArticleGoogle Scholar
- Austin DW, et al., et al.: Gene network shaping of inherent noise spectra. Nature 2006,439(7076):608-611. 10.1038/nature04194View ArticleGoogle Scholar
- Becskei A, Serrano L: Engineering stability in gene networks by autoregulation. Nature 2000,405(6786):590-593. 10.1038/35014651View ArticleGoogle Scholar
- Savageau M: Comparison of classical and autogenous systems of regulation in inducible operons. Nature 1974, 252: 546-549. 10.1038/252546a0View ArticleGoogle Scholar
- Zhang Z, Qian W, Zhang J: Positive selection for elevated gene expression noise in yeast. Mol Syst Biol 2009, 5: 299.View ArticleGoogle Scholar
- Bar-Even A, et al., et al.: Noise in protein expression scales with natural protein abundance. Nat Genet 2006,38(6):636-643. 10.1038/ng1807View ArticleGoogle Scholar
- Newman JRS, et al., et al.: Single-cell proteomic analysis of S. cerevisiae reveals the architecture of biological noise. Nature 2006,441(7095):840-846. 10.1038/nature04785View ArticleGoogle Scholar
- Fraser H, et al., et al.: Noise minimization in eukaryotic gene expression. PLoS Biol 2004,2(6):834-838.View ArticleGoogle Scholar
- Dublanche Y, et al., et al.: Noise in transcription negative feedback loops: simulation and experimental analysis. Mol Syst Biol 2006, 2: 41.View ArticleGoogle Scholar
- Nevozhay D, et al., et al.: Negative autoregulation linearizes the dose–response and suppresses the heterogeneity of gene expression. Proc Natl Acad Sci U S A 2009,106(13):5123-5128. 10.1073/pnas.0809901106View ArticleGoogle Scholar
- Simpson ML, Cox CD, Sayler GS: Frequency domain analysis of noise in autoregulated gene circuits. Proc Natl Acad Sci U S A 2003,100(8):4551-4556. 10.1073/pnas.0736140100View ArticleGoogle Scholar
- Thattai M, Shraiman BI: Metabolic switching in the sugar phosphotransferase system of Escherichia coli. Biophys J 2003,85(2):744-754. 10.1016/S0006-3495(03)74517-2View ArticleGoogle Scholar
- Wolf DM, Vazirani VV, Arkin AP: Diversity in times of adversity: probabilistic strategies in microbial survival games. J Theor Biol 2005,234(2):227-253. 10.1016/j.jtbi.2004.11.020View ArticleGoogle Scholar
- Brock A, Chang H, Huang S: Non-genetic heterogeneity–a mutation-independent driving force for the somatic evolution of tumours. Nat Rev Genet 2009,10(5):336-342. 10.1038/nrg2556View ArticleGoogle Scholar
- Charlebois DA, Abdennur N, Kaern M: Gene expression noise facilitates adaptation and drug resistance independently of mutation. Phys Rev Lett 2011,107(21):218101.View ArticleGoogle Scholar
- Fraser D, Kaern M: A chance at survival: gene expression noise and phenotypic diversification strategies. Mol Microbiol 2009,71(6):1333-1340. 10.1111/j.1365-2958.2009.06605.xView ArticleGoogle Scholar
- Zhuravel D, et al., et al.: Phenotypic impact of regulatory noise in cellular stress-response pathways. Syst Synth Biol 2010,4(2):105-116. 10.1007/s11693-010-9055-2View ArticleGoogle Scholar
- Lewis K: Persister cells, dormancy and infectious disease. Nat Rev Microbiol 2007,5(1):48-56. 10.1038/nrmicro1557View ArticleGoogle Scholar
- Gerdes K, Christensen S, Lobner-Olesen A: Prokaryotic toxin-antitoxin stress response loci. Nat Rev Microbiol 2005,3(5):371-382. 10.1038/nrmicro1147View ArticleGoogle Scholar
- Lewis K: Persister cells. Annu Rev Microbiol 2010, 64: 357-372. 10.1146/annurev.micro.112408.134306View ArticleGoogle Scholar
- Bigger J: Treatment of staphylococcal infections with penicillin by intermittent sterilisation. Lancet 1944, ii: 497-500.View ArticleGoogle Scholar
- Christensen S, Gerdes K: RelE toxins from Bacteria and Archaea cleave mRNAs on translating ribosomes, which are rescued by tmRNA. Mol Microbiol 2003,48(5):1389-1400. 10.1046/j.1365-2958.2003.03512.xView ArticleGoogle Scholar
- Lou C, Li Z, Ouyang Q: A molecular model for persister in E. coli. J Theor Biol 2008,255(2):205-209. 10.1016/j.jtbi.2008.07.035View ArticleGoogle Scholar
- Ramage HR, Connolly LE, Cox JS: Comprehensive Functional Analysis of Mycobacterium tuberculosis Toxin-Antitoxin Systems: Implications for Pathogenesis, Stress Responses, and Evolution. PLoS Genet 2009,5(12):e1000767. 10.1371/journal.pgen.1000767View ArticleGoogle Scholar
- Schumacher MA, et al., et al.: Molecular Mechanisms of HipA-Mediated Multidrug Tolerance and Its Neutralization by HipB. Science 2009,323(5912):396-401. 10.1126/science.1163806View ArticleGoogle Scholar
- Balaban N, et al., et al.: Bacterial persistence as a phenotypic switch. Science 2004,305(5690):1622-1625. 10.1126/science.1099390View ArticleGoogle Scholar
- Black D, Irwin B, Moyed H: Autoregulation of Hip, an operon that affects lethality due to inhibition of peptidoglycan of DNA-synthesis. J Bacteriol 1994,176(13):4081-4091.Google Scholar
- Vilar JMG, et al., et al.: Mechanisms of noise-resistance in genetic oscillators. Proc Natl Acad Sci U S A 2002,99(9):5988-5992. 10.1073/pnas.092133899View ArticleGoogle Scholar
- Rotem E, et al., et al.: Regulation of phenotypic variability by a threshold-based mechanism underlies bacterial persistence. Proc Natl Acad Sci U S A 2010,107(28):12541-12546. 10.1073/pnas.1004333107View ArticleGoogle Scholar
- Bernstein J, et al., et al.: Global analysis of Escherichia coli RNA degradosome function using DNA microarrays. Proc Natl Acad Sci U S A 2004,101(9):2758-2763. 10.1073/pnas.0308747101View ArticleGoogle Scholar
- Gillespie DT: Exact Stochastic Simulation Of Coupled Chemical-Reactions. J Phys Chem 1977,81(25):2340-2361. 10.1021/j100540a008View ArticleGoogle Scholar
- Acar M, Becskei A, van Oudenaarden A: Enhancement of cellular memory by reducing stochastic transitions. Nature 2005,435(7039):228-232. 10.1038/nature03524View ArticleGoogle Scholar
- Choi PJ, et al., et al.: A stochastic single-molecule event triggers phenotype switching of a bacterial cell. Science 2008,322(5900):442-446. 10.1126/science.1161427View ArticleGoogle Scholar
- Dhar N, McKinney JD: Microbial phenotypic heterogeneity and antibiotic tolerance. Curr Opin Microbiol 2007,10(1):30-38. 10.1016/j.mib.2006.12.007View ArticleGoogle Scholar
- Kussell E, et al., et al.: Bacterial persistence: A model of survival in changing environments. Genetics 2005,169(4):1807-1814. 10.1534/genetics.104.035352View ArticleGoogle Scholar
- Murray JD: Mathematical Biology I. Springer, ; 2002.Google Scholar