Modeling lymphocyte homing and encounters in lymph nodes
© Baldazzi et al; licensee BioMed Central Ltd. 2009
Received: 10 March 2009
Accepted: 25 November 2009
Published: 25 November 2009
The efficiency of lymph nodes depends on tissue structure and organization, which allow the coordination of lymphocyte traffic. Despite their essential role, our understanding of lymph node specific mechanisms is still incomplete and currently a topic of intense research.
In this paper, we present a hybrid discrete/continuous model of the lymph node, accounting for differences in cell velocity and chemotactic response, influenced by the spatial compartmentalization of the lymph node and the regulation of cells migration, encounter, and antigen presentation during the inflammation process.
Our model reproduces the correct timing of an immune response, including the observed time delay between duplication of T helper cells and duplication of B cells in response to antigen exposure. Furthermore, we investigate the consequences of the absence of dendritic cells at different times during infection, and the dependence of system dynamics on the regulation of lymphocyte exit from lymph nodes. In both cases, the model predicts the emergence of an impaired immune response, i.e., the response is significantly reduced in magnitude. Dendritic cell removal is also shown to delay the response time with respect to normal conditions.
Lymph nodes and Peyer's patches play key roles in the development of an appropriate and efficient immune response. Once an Antigen (Ag) is captured by Ag-processing cells, it is rapidly carried to the nearest lymph node, where it is presented to specific lymphocytes to trigger an immune response. The recognition phase must be highly efficient: within a few hours, it is necessary to find specific lymphocytes among a repertoire that includes a very large number of receptors [1, 2]. The specific architecture of the lymph node and a fine-tuned balance between diffusion, chemotaxis, and receptor expression are the basis of this process.
Human lymph nodes are bean-shaped structures that range in size from a few millimeters to about 1-2 cm in their normal state. Internally, two main regions can be distinguished: the medulla and the cortex. The cortex can be further divided into an inner part, the paracortex (also called the T cell area), rich in T lymphocytes and an outer area, the node cortex that includes the B cell area consisting of follicles and germinal centers, where B cells are activated and differentiate . T and B areas are identified by high concentrations of different chemokines (CCR7 and CXCR5, respectively) secreted by local stromal cells [1, 4, 5]. The whole structure is supported by a dense network of fibroblastic reticular cells that encloses small lymphatic channels of 10-15 μ m in diameter along which small molecules are thought to diffuse . Macrophages, dendritic cells (DC), and some lymphocytes flow from the afferent lymphatic vessels, through the fibroblastic reticular cellular network, to the node cortex and the medulla, before leaving via the efferent lymphatic vessels. Most T and B cells, however, are too large in size and enter the lymph node mainly from the blood, through high endothelial venules (HEV) located inside the paracortex. During an infection, lymphocyte recruitment from the periphery is enhanced due to a widening of the primary arteriole feeding the lymph node . Once inside the lymph node, B and T cells rapidly home in their own compartments, following a specific chemotactic gradient [1, 4]. T helper cells (TH) are the fastest, with an average velocity of 11 μ m/min, followed by B cells with 6 μ m/min and DCs with a velocity of 3 μ m/min .
In the absence of an antigenic challenge, T and B cells randomly scan their respective areas for ~24-48 h before exiting the lymph node . Entrance of antigens into the lymph node triggers a series of events leading to antigen recognition and the activation of an immune response.
Two distinct pathways for antigen delivery have been recognized. In general, an Ag-presenting cell picks up the antigen in peripheral tissues and migrates to the nearest lymph node in order to present the MHC-peptide complex to the T and B cells. In addition, antigens of low molecular weight (below 70 kDa) may enter the lymph node as soluble antigens and reach the T cell area directly through the fibroblastic reticular cellular meshwork, without any preliminary recognition. This delivery system is very efficient, and soluble antigens are detected inside the lymph node within a few minutes of infection. Here, specific resident DCs are able to capture the antigen well before Ag-presenting cells arrive from peripheral organs (on average, 8-12 h later) .
The immune response is initiated inside the T cell area, where antigens are first presented to the T cells, generally by DCs. After antigen recognition, peripheral dendritic cells undergo a change in the expression of their surface receptors. In particular, receptors for inflammatory chemokines are lost, and lymphoid receptors, especially the CCR7 receptor, are expressed . As a consequence, Ag-presenting DCs are rapidly routed towards the T cell area, where the probability of meeting specific T cells is higher. In addition, once inside the lymph node, Ag-loaded DCs begin to release specific chemokines (e.g., MDC) to guide proximal active T cells [11, 12].
B cell activation is initiated by an engagement of the B cell receptor either by a soluble or a membrane-associated antigen [13, 14]. Following antigenic stimulation, B cells start to co-express the CCR7 receptor [15, 16] and rapidly localize at the boundary between the T and B areas. Here, they can receive the right costimulation by T helper cells and start their proliferation and differentiation processes. Unstimulated lymphocytes rapidly pass through the lymph node to return to general circulation.
Lymphocyte egress from lymph nodes is still a subject of investigation. Recent studies suggest a key role of the Sphingosine-1-phosphate (S1P) molecule [17–19]. Its specific receptor S1P1 acts as a type of "pass filter" that selectively controls lymphocyte exit through efferent lymphatic vessels . The receptor S1P1 is downregulated during lymphocyte activation , preventing Ag-specific lymphocytes from leaving the lymph node before an immune response has been mounted.
The variety of mechanisms and the large number of entities present in a lymph node make it difficult to understand the role played by each single component in the overall behavior. Models for several different aspects of lymphocyte motion have been developed [20–22]. Here, we present a model of the lymph node that is able to capture the interplay between different mechanisms, and which assures a fast and efficient immune response. We focus on lymphocyte recruitment and trafficking inside the lymph node, including specific cell diffusion properties, chemotaxis, and control of lymphocyte egress. We resort to a hybrid discrete/continuous approach that combines a stochastic agent-based description of cell interactions with a continuous model of molecular diffusion described by partial differential equations.
Agent-based models have a long history in immunology [23–25] because they have proven to be well-suited to handle complex nonlinearities and differences in individual cell characteristics. The choice of a hybrid approach is motivated by the need to explicitly consider the effects of chemotaxis on cells motility. This requires the introduction of a new timescale and, as a consequence, a distinct level of representation with respect to the inter-cellular interactions.
The model currently includes the description of three types of cells: TH cells, B cells, and dendritic cells, plus four chemokines, CXCR5, CCR7, S1P, and MDC.
Simple molecules like chemokines are represented by their spatial concentration that acts as a signal that triggers cellular interactions and movement. Unlike cells, molecules are not endowed with an internal description.
Cellular internal states.
Cells change their internal state when they interact with other entities. Interactions are coded as probabilistic rules; each interaction requires the cells involved to be in a specific state and the interaction probability depends on the complementarity of their receptor strings, as measured by their Hamming distance. The interactions are local and many entities may sit on each single lattice site. In particular, for each lattice point, all interaction rules are executed in random order and a greedy paradigm is applied (i.e., cells that have interacted are taken out from the pool for the current time step). The list of all included interactions and signaling mechanisms is shown in Fig. 2.
The model description of chemotaxis is difficult due to the range of timescales involved in the phenomenon. Cells, in fact, are characterized by a diffusion coefficient much lower than for small signaling molecules , that can rapidly diffuse throughout the entire simulation volume.
where c = c(x) is the concentration of chemokines (expressed in pM), D is the diffusion coefficient, and λ is proportional to the half-life of the molecules. We assume that D = 3000 μ m2/min and λ = 3 h [25, 26]. Homing chemokines are continuously released by several sources that are randomly distributed inside the T and B areas. At each time step, a burst of new molecules is injected at each source point and immediately spread by diffusion over the simulation volume.
In contrast to chemokines, cells move individually at discrete time steps (one simulation time step ~30 mins real time). Chemotaxis is modeled in the following way. First, we compute the normalized difference of chemoattractant concentration between the current cell position and the neighboring lattice points, yielding a set of discrete percentages of total neighboring chemoattractant. The probability of moving in each direction is then partitioned according to the resulting percentages, and the direction of cell displacement is chosen stochastically. This results in cells moving with a higher frequency towards the maximum increase of chemoattractant concentration, and with a lower frequency away from chemoattractants, although this possibility is allowed. However, a minimal intensity for the chemotactic gradient is required in order to elicit a chemotactic response; a weak or null chemotactic field results in a simple random walk, similar to the one observed in two-photon microscopy imaging experiments [27, 28]. The threshold is assumed to be cell-dependent. T cells are provided with a lower chemotactic sensitivity than B and DC cells (of a factor 103) . This allows T cells to move essentially randomly inside their area whereas B and dendritic cells are progressively guided towards the appropriate position (respectively, the B cell area and the T cell area) during the immune response [3, 29].
Differences in cell mobility are also taken into account. Within the limit of our temporal resolution, the dynamics of cell-cell interactions is not relevant, and cell velocity is assumed to be constant during the entire simulation.
The transient effects of the chemotactic signals on a given cell are modeled by a modulation in the expression of the receptor signal, according to the cell's internal state. For instance, when a B cell presents the antigen, the expression of its CXCR5 receptor is turned off, and the CCR7 receptor is expressed in its place. In this way, B cells are redirected towards the boundary between the T and B cell areas, to receive the appropriate costimulation.
where D TH is the diffusion coefficient for the TH cell and Δt is the time step.
At each time step, a constant number of naïve TH and B cells enter on the average the lymph node through several high endothelial venules, randomly distributed inside the T cell area . Dendritic cells enter, instead, through the lymphatic flux, from the afferent vessel. A continuous renewal of lymphocytes and dendritic cells is provided by assuming a total population of about 106 cells (fluxes have been taken from ). To model the increase in lymphocyte recruitment during an infection, the incoming flux of TH and B cells is augmented by a factor that is proportional to the number of antigens inside the lymph node. We represent antigens as small immunogenic molecules. The model includes the description of both mechanisms for antigen delivery, based on a passive or active (i.e., DC-aided) transport mechanism. Soluble Ags are introduced into the lymph node immediately after infection, at random sites inside the T compartment. DC-presenting cells, instead, are inserted in the simulation volume shortly after the injection of soluble Ags, approximately 8-12 h later.
Parameter set - parameters of the simulation.
B's duplication steps
TH's duplication steps
B cells' velocity
6 μ m/min
11 μ m/min
3 μ m/min
chemokines' diffusion coeff.
3·103 μ m2/min
chemokines' half life
B's initials per μ l
105 (cells/μ l)
TH's initials per μ l
2·105 (cells/μ l)
DC's initials per μ l
103 (cells/μ l)
incoming flux of TH cells
incoming flux of TH cells
DC affinity to Ag (prob)
numbers of HEV per μ l
4 (μ l-1)
exit prob for cells
released chemokines qty
108 mol (min src)-1
10 (μ l)
time steps of the simulation
injection freq. of soluble Ags
soluble Ags per injection
The tuning of free parameters was performed by comparing the results against experimental values taken from the literature. For instance, i) an average lymphocyte residence time inside the lymph node of approximately 24-48 h ; ii) an outgoing flux per day of 25% of the population of cells ; (iii) an average homing time for newly-arrived B lymphocytes of about 10 hrs (estimated from ).
Results and Discussion
Lymph node function is characterized as being in one of two distinct regimes or states, namely healthy or infected. Each state, and the transition between states, is characterized by specific timescales that result from the delicate combination of structural architecture, diffusion properties of components, and flux remodeling, as described above. Our model succeeds in reproducing the expected dynamics of the system with and without the presence of antigenic stimulation.
As long as new antigens are detected, the immune response persists: the total lymphocyte number increases monotonically up to 5 times the original value [2, 8], and cells continue to duplicate, increasing the pool of available Ag-specific cells (see Fig. 5). Once the infection is over, the immune system and the lymph node recover the original healthy state. After about 200 h, in the presented simulation, Ag concentration has nearly vanished, and the incoming cell flux begins a decrease to its original value. The total number of TH and B cells gradually reduces and recovers its normal value within a couple of weeks. B cells leave the lymph node at a lower rate because of their reduced motility with respect to TH cells.
In silico experiments
We now use the model to investigate different scenarios by perturbing the normal function of the system and looking at the emergent response. Several experiments can be designed in response to specific questions. In the following, we propose two distinct experiments aimed at investigating the biological role of DC and S1P1-control of lymphocyte egress on the onset of an efficient immune response.
The role of dendritic cells
Dendritic cells are recognized as the most efficient antigen-presenting cells. Their function is of primary importance for the activation of specific lymphocytes inside a lymph node. In the following, we show the effect of removing all DCs from the lymph node at a given time. Depending on the delay from Ag injection, the effect on the immune response varies.
Lymphocyte exit control: role of the S1P receptor
The mechanisms that regulate lymphocyte exit from the lymph node have been subject to intense research activity in the last few years. Several studies implicate the Sphingosine-1-phosphate receptor (S1P1) as a crucial element for the regulation of the exit mechanism of lymphocytes from the lymph node and highlight its importance in the selective retention of Ag-specific lymphocytes during an infection.
However, the S1P1 control does not seem to affect the onset and overall dynamics of the immune response: the timing of lymphocyte encounters appears to be mainly determined by the geometry and the diffusion properties of the system rather than the Ag-specificity of its cell population.
We have developed a model that reproduces some aspects of the immune response and the behavior of cell/antigen motility within a lymph node. In particular, we focus on the mechanisms that determine the onset of a primary immune response, from Ag delivery to B cell activation and duplication and the observed lymph node shrinking after an immune response. We stress that the obtained T and B cell responses at the correct times are emergent properties of a quasi-realistic description of lymphocyte density, interactions, and motion (including chemotaxis and diffusion characteristics), combined with a schematic description of lymph node compartmentalization. This is a key difference from previous work  in which the correct timing of the immune response was somehow hard-coded in the simulator.
The model provides interesting insights into the role played by DCs and by the regulation of lymphocyte exit from lymph nodes on the resulting immune response. We show an impaired immune response when one of these mechanisms is perturbed. DC removal at early times produces large effects, with an immune response that is greatly delayed and reduced in magnitude due to the lack of active TH cells that are able to provide the right costimulation of B lymphocytes. S1P1 control of specific cells, instead, affects essentially the magnitude of the immune response rather than the timing, by decreasing the overall number of specific lymphocytes that can participate in Ag detection. In the same spirit, many other experiments could be planned to investigate, for example, the role of chemotaxis and transient receptor expression on the different phases of Ag presentation.
The current model represents a first attempt to comprehensively sketch fundamental aspects of lymph node function based on a few cellular component and molecular mechanisms. Only essential pathways have been included in the current version of the model, and further work is needed to enhance the description of the rich variety of behaviors that can be observed in a real immune response. Future improvements include a better description of B cell activation (by macrophages or dendritic cells) and differentiation into plasma and memory cells. In particular, it would be interesting to examine the impact of T-independent B activation mechanisms on the emergence of a humoral immune response . Moreover, the presence of memory lymphocytes has important consequences for the dynamics of a second immunization. With the inclusion of memory cell generation in the model, differences in lymphocyte population and Ag-presentation efficiency between the primary and secondary responses could be analyzed. Specific Ag features (e.g., if it is a bacterium or a virus) and other related biological processes could also be taken into account to address specific questions. In parallel with progress in lymph node understanding, this simulator can be a useful tool to test new hypotheses, investigating the effect of additional mechanisms on the resulting immune response.
Availability and requirements
We thank the "Consorzio interuniversitario per le Applicazioni di Supercalcolo Per Università e Ricerca" (CASPUR) for computing resources and support. We also thank Mark Coles for useful and clarifying discussions. This work was partially supported by the European Community under the EC contract FP6-2004-IST-4, No.028069 (ImmunoGrid).
- Cyster J: Chemokines and cell migration in secondary lymphoid organs. Science 1999, 286: 2098–2102. 10.1126/science.286.5447.2098View ArticlePubMedGoogle Scholar
- Soderberg KA, Payne GW, Sato A, Medzhitov R, Segal SS, Iwasaki A: Innate control of adaptive immunity via remodeling of lymph node feed arteriole. Proc Natl Acad Sci USA 2005, 102: 16315–16320. 10.1073/pnas.0506190102PubMed CentralView ArticlePubMedGoogle Scholar
- Garside P, Ingulli E, Merica RR, Johnson JG, Noelle RJ, Jenkins MK: Visualization of Specific B and T Lymphocyte Interactions in the Lymph Node. Science 1998, 281(5373):96–99. 10.1126/science.281.5373.96View ArticlePubMedGoogle Scholar
- Moser B, Loetscher P: Lymphocyte traffic control by chemokines. Nat Immunol 2001, 2: 123–128. 10.1038/84219View ArticlePubMedGoogle Scholar
- Catron DM, Itano AA, Pape KA, Mueller DL, Jenkins MK: Visualizing the First 50 Hr of the Primary Immune Response to a Soluble Antigen. Immunity 2004, 21: 341–347. 10.1016/j.immuni.2004.08.007View ArticlePubMedGoogle Scholar
- Willard-Mack CL: Normal structure, function, and histology of lymph nodes. Toxicol Pathol 2006, 34(5):409–24. 10.1080/01926230600867727View ArticlePubMedGoogle Scholar
- Miller MJ, Wei SH, Cahalan MD, Parker I: T cell repertoire scanning is promoted by dynamic dendritic cell behavior and random T cell motility in the lymph node. Proc Natl Acad Sci USA 2004, 101: 998–1003. 10.1073/pnas.0306407101PubMed CentralView ArticlePubMedGoogle Scholar
- Qi H, Egen JG, Huang AYC, Germain RN: Extrafollicular Activation of Lymph Node B Cells by Antigen-Bearing Dendritic Cells. Science 2006, 312(5780):1672–1676. 10.1126/science.1125703View ArticlePubMedGoogle Scholar
- Sixt M, Kanazawa N, Selg M, Samson T, Roos G, Reinhardt DP, Pabst R, Lutz MB, Sorokin L: The Conduit System Transports Soluble Antigens from the Afferent Lymph to Resident Dendritic Cells in the T Cell Area of the Lymph Node. Immunity 2005, 22: 19–29. 10.1016/j.immuni.2004.11.013View ArticlePubMedGoogle Scholar
- von Andrian UH, Mempel TR: Homing and cellular traffic in lymph nodes. Nat Rev Immunol 2003, 3(11):867–878. 10.1038/nri1222View ArticlePubMedGoogle Scholar
- Katou F, Ohtani H, Nakayama T, Ono K, Matsushima K, Saaristo A, Nagura H, Yoshie O, Motegi K: Macrophage-Derived Chemokine (MDC/CCL22) and CCR4 Are Involved in the Formation of T Lymphocyte-Dendritic Cell Clusters in Human Inflamed Skin and Secondary Lymphoid Tissue. Am J Pathol 2001, 158(4):1263–1270.PubMed CentralView ArticlePubMedGoogle Scholar
- Kaiser A, Donnadieu E, Abastado JP, Trautmann A, Nardin A: CC Chemokine Ligand 19 Secreted by Mature Dendritic Cells Increases Naive T Cell Scanning Behavior and Their Response to Rare Cognate Antigen. J Immunol 2005, 175(4):2349–2356.View ArticlePubMedGoogle Scholar
- Bergtold A, Desai DD, Gavhane A, Clynes R: Cell Surface Recycling of Internalized Antigen Permits Dendritic Cell Priming of B Cells. Immunity 2005, 23(5):503–514. 10.1016/j.immuni.2005.09.013View ArticlePubMedGoogle Scholar
- Okada T, Cyster J: B cell migration and interactions in the early phase of antibody responses. Curr Opin Immunol 2006, 18(3):278–285. 10.1016/j.coi.2006.02.005View ArticlePubMedGoogle Scholar
- Reif K, Ekland E, Ohl L, Nakano H, Lipp M, Förster R, Cyster J: Balanced responsiveness to chemoattractants from adjacent zones determines B-cell position. Nature 2002, 416: 94–99. 10.1038/416094aView ArticlePubMedGoogle Scholar
- Hardtke S, Ohl L, Forster R: Balanced expression of CXCR5 and CCR7 on follicular T helper cells determines their transient positioning to lymph node follicles and is essential for efficient B-cell help. Blood 2005, 106(6):1924–1931. 10.1182/blood-2004-11-4494View ArticlePubMedGoogle Scholar
- Matloubian M, Lo CG, Cinamon G, Lesneski MJ, Xu Y, Brinkmann V, Allende ML, Proia RL, Cyster JG: Lymphocyte egress from thymus and peripheral lymphoid organs is dependent on S1P receptor 1. Nature 2004, 427(6972):355–360. 10.1038/nature02284View ArticlePubMedGoogle Scholar
- Allende ML, Dreier JL, Mandala S, Proia RL: Expression of the Sphingosine 1-Phosphate Receptor, S1P1, on T-cells Controls Thymic Emigration. J Biol Chem 2004, 279(15):15396–15401. 10.1074/jbc.M314291200View ArticlePubMedGoogle Scholar
- Wei SH, Rosen H, Matheu MP, Sanna GM, Wang SK, Jo E, Wong CH, Parker I, Cahalan MD: Sphingosine 1-phosphate type 1 receptor agonism inhibits transendothelial migration of medullary T cells to lymphatic sinuses. Nat Immunol 2005, 6(12):1228–1235. 10.1038/ni1269View ArticlePubMedGoogle Scholar
- Kesmir C, De Boer RJ: A Mathematical Model on Germinal Center Kinetics and Termination. J Immunol 1999, 163(5):2463–2469.PubMedGoogle Scholar
- Kleinstein SH, Singh JP: Toward Quantitative Simulation of Germinal Center Dynamics: Biological and Modeling Insights from Experimental Validation. J Theor Biol 2001, 211: 253–275. 10.1006/jtbi.2001.2344View ArticlePubMedGoogle Scholar
- Beltman JB, Maree AF, de Boer RJ: Spatial modelling of brief and long interactions between T cells and dendritic cells. Immunol Cell Biol 2007, 85: 306–314. 10.1038/sj.icb.7100054View ArticlePubMedGoogle Scholar
- Celada F, Seiden P: A computer model of cellular interaction in the immune system. Immunol Today 1992, 13: 56–62. 10.1016/0167-5699(92)90135-TView ArticlePubMedGoogle Scholar
- Bernaschi M, Castiglione F: Design and implementation of an immune system simulator. Comput Biol Med 2001, 31: 303–331. 10.1016/S0010-4825(01)00011-7View ArticlePubMedGoogle Scholar
- Segovia-Juarez JL, Ganguli S, Kirschner D: Identifying control mechanisms of granuloma formation during M. tuberculosis infection using an agent-based model. J Theor Biol 2004, 231(3):357–376. 10.1016/j.jtbi.2004.06.031View ArticlePubMedGoogle Scholar
- Francis K, Palsson BO: Effective intercellular communication distances are determined by the relative time constants for cyto/chemokine secretion and diffusion. Proc Natl Acad Sci USA 1997, 94(23):12258–12262. 10.1073/pnas.94.23.12258PubMed CentralView ArticlePubMedGoogle Scholar
- Miller MJ, Wei SH, Parker I, Cahalan MD: Two-Photon Imaging of Lymphocyte Motility and Antigen Response in Intact Lymph Node. Science 2002, 296(5574):1869–1873. 10.1126/science.1070051View ArticlePubMedGoogle Scholar
- Miller MJ, Wei SH, Cahalan MD, Parker I: Autonomous T cell trafficking examined in vivo with intravital two-photon microscopy. Proc Natl Acad Sci USA 2003, 100: 2604–2609. 10.1073/pnas.2628040100PubMed CentralView ArticlePubMedGoogle Scholar
- Okada T, Miller M, Parker I, Krummel M, Neighbors Ma: Antigen-Engaged B Cells Undergo Chemotaxis toward the T Zone and Form Motile Conjugates with Helper T Cells. Plos Biol 2005, 3(6):1047–1061. 10.1371/journal.pbio.0030150View ArticleGoogle Scholar
- Wei SH, Parker I, Miller MJ, Cahalan MD: A stochastic view of lymphocyte motility and trafficking within the lymph node. Immunol Rev 2003, 195: 136–159. 10.1034/j.1600-065X.2003.00076.xView ArticlePubMedGoogle Scholar
- Young AJ: The physiology of lymphocyte migration through the single lymph node in vivo. Seminars in Immunology 1999, 11(2):73–83. 10.1006/smim.1999.0163View ArticlePubMedGoogle Scholar
- Mempel TR, Henrickson SE, Von Andrian UH: T-cell priming by dendritic cells in lymph nodes occurs in three distinct phases. Nature 2004, 427(6970):154–159. 10.1038/nature02238View ArticlePubMedGoogle Scholar
- Fagarasan S, Honjo T: T-Independent Immune Response: New Aspects of B Cell Biology. Science 2000, 290(5489):89–92. 10.1126/science.290.5489.89View ArticlePubMedGoogle Scholar
- Goldsby R, TJ TK, Osborne B: Kuby Immunology. IV edition. New York: W.H. Freeman and Company; 2003.Google Scholar
- Murray JD: Mathematical Biology. Volume II. springer; 2003.Google Scholar
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.