Applications of mechanistic modelling to clinical and experimental immunology: an emerging technology to accelerate immunotherapeutic discovery and development.

The application of in-silico modelling is beginning to emerge as a key methodology to advance our understanding of mechanisms of disease pathophysiology and related drug action, and in the design of experimental medicine and clinical studies. From this perspective, we will present a non-technical discussion of a small number of recent and historical applications of mathematical, statistical and computational modelling to clinical and experimental immunology. We focus specifically upon mechanistic questions relating to human viral infection, tumour growth and metastasis and T cell activation. These exemplar applications highlight the potential of this approach to impact upon human immunology informed by ever-expanding experimental, clinical and 'omics' data. Despite the capacity of mechanistic modelling to accelerate therapeutic discovery and development and to de-risk clinical trial design, it is not utilized widely across the field. We outline ongoing challenges facing the integration of mechanistic modelling with experimental and clinical immunology, and suggest how these may be overcome. Advances in key technologies, including multi-scale modelling, machine learning and the wealth of 'omics' data sets, coupled with advancements in computational capacity, are providing the basis for mechanistic modelling to impact on immunotherapeutic discovery and development during the next decade.

Mechanical Cell-Cell Communication in Fibrous Networks: The Importance of Network Geometry.

Cells contracting in extracellular matrix (ECM) can transmit stress over long distances, communicating their position and orientation to cells many tens of micrometres away. Such phenomena are not observed when cells are seeded on substrates with linear elastic properties, such as polyacrylamide (PA) gel. The ability for fibrous substrates to support far reaching stress and strain fields has implications for many physiological processes, while the mechanical properties of ECM are central to several pathological processes, including tumour invasion and fibrosis. Theoretical models have investigated the properties of ECM in a variety of network geometries. However, the effects of network architecture on mechanical cell-cell communication have received little attention. This work investigates the effects of geometry on network mechanics, and thus the ability for cells to communicate mechanically through different networks. Cell-derived displacement fields are quantified for various network geometries while controlling for network topology, cross-link density and micromechanical properties. We find that the heterogeneity of response, fibre alignment, and substrate displacement fields are sensitive to network choice. Further, we show that certain geometries support mechanical communication over longer distances than others. As such, we predict that the choice of network geometry is important in fundamental modelling of cell-cell interactions in fibrous substrates, as well as in experimental settings, where mechanical signalling at the cellular scale plays an important role. This work thus informs the construction of theoretical models for substrate mechanics and experimental explorations of mechanical cell-cell communication.

Vascular phenotype identification and anti-angiogenic treatment recommendation: A pseudo-multiscale mathematical model of angiogenesis.

The development of anti-angiogenic drugs for cancer therapy has yielded some promising candidates, but novel approaches for interventions to angiogenesis have led to disappointing results. In addition, there is a shortage of biomarkers that are predictive of response to anti-angiogenic treatments. Consequently, the complex biochemical and physiological basis for tumour angiogenesis remains incompletely understood. We have adopted a mathematical approach to address these issues, formulating a spatially averaged multiscale model that couples the dynamics of VEGF, Ang1, Ang2 and PDGF, with those of mature and immature endothelial cells and pericyte cells. The model reproduces qualitative experimental results regarding pericyte coverage of vessels after treatment by anti-Ang2, anti-VEGF and combination anti-VEGF/anti-Ang2 antibodies. We used the steady state behaviours of the model to characterise angiogenic and non-angiogenic vascular phenotypes, and used mechanistic perturbations representing hypothetical anti-angiogenic treatments to generate testable hypotheses regarding transitions to non-angiogenic phenotypes that depend on the pre-treatment vascular phenotype. Additionally, we predicted a synergistic effect between anti-VEGF and anti-Ang2 treatments when applied to an immature pre-treatment vascular phenotype, but not when applied to a normalised angiogenic pre-treatment phenotype. Based on these findings, we conclude that changes in vascular phenotype are predicted to be useful as an experimental biomarker of response to treatment. Further, our analysis illustrates the potential value of non-spatial mathematical models for generating tractable predictions regarding the action of anti-angiogenic therapies.

Investigating the Turing conditions for diffusion-driven instability in the presence of a binding immobile substrate.

Turing's diffusion-driven instability for the standard two species reaction-diffusion system is only achievable under well-known and rather restrictive conditions on both the diffusion rates and the kinetic parameters, which necessitates the pairing of a self-activator with a self-inhibitor. In this study we generalize the standard two-species model by considering the case where the reactants can bind to an immobile substrate, for instance extra-cellular matrix, and investigate the influence of this dynamics on Turing's diffusion-driven instability. Such systems have been previously studied on the grounds that binding of the self-activator to a substrate may effectively reduce its diffusion rate and thus induce a Turing instability for species with equal diffusion coefficients, as originally demonstrated by Lengyel and Epstein (1992) under the assumption that the bound state dynamics occurs on a fast timescale. We, however, analyse the full system without any separation of timescales and demonstrate that the full system also allows a relaxation of the standard constraints on the reaction kinetics for the Turing instability, increasing the type of interactions that could give rise to spatial patterning. In particular, we show that two self-activators can undertake a diffusively driven instability in the presence of a binding immobile substrate, highlighting that the interactions required of a putative biological Turing instability need not be associated with a self-activator-self-inhibitor morphogen pair.

Modelling Aedes aegypti mosquito control via transgenic and sterile insect techniques: endemics and emerging outbreaks.

The invasion of pest insects often changes or destroys a native ecosystem, and can result in food shortages and disease endemics. Issues such as the environmental effects of chemical control methods, the economic burden of maintaining control strategies and the risk of pest resistance still remain, and mosquito-borne diseases such as malaria and dengue fever prevail in many countries, infecting over 100 million worldwide in 2010. One environmentally friendly method for mosquito control is the Sterile Insect Technique (SIT). This species-specific method of insect control relies on the mass rearing, sterilization and release of large numbers of sterile insects. An alternative transgenic method is the Release of Insects carrying a Dominant Lethal (RIDL). Our objective is to consider contrasting control strategies for two invasive scenarios via SIT and RIDL: an endemic case and an emerging outbreak. We investigate how the release rate and size of release region influence both the potential for control success and the resources needed to achieve it, under a range of conditions and control strategies, and we discuss advantageous strategies with respect to reducing the release resources and strategy costs (in terms of control mosquito numbers) required to achieve complete eradication of wild-type mosquitoes.

Modelling biological invasions: Individual to population scales at interfaces.

Extracting the population level behaviour of biological systems from that of the individual is critical in understanding dynamics across multiple scales and thus has been the subject of numerous investigations. Here, the influence of spatial heterogeneity in such contexts is explored for interfaces with a separation of the length scales characterising the individual and the interface, a situation that can arise in applications involving cellular modelling. As an illustrative example, we consider cell movement between white and grey matter in the brain which may be relevant in considering the invasive dynamics of glioma. We show that while one can safely neglect intrinsic noise, at least when considering glioma cell invasion, profound differences in population behaviours emerge in the presence of interfaces with only subtle alterations in the dynamics at the individual level. Transport driven by local cell sensing generates predictions of cell accumulations along interfaces where cell motility changes. This behaviour is not predicted with the commonly used Fickian diffusion transport model, but can be extracted from preliminary observations of specific cell lines in recent, novel, cryo-imaging. Consequently, these findings suggest a need to consider the impact of individual behaviour, spatial heterogeneity and especially interfaces in experimental and modelling frameworks of cellular dynamics, for instance in the characterisation of glioma cell motility.

Meniscal tear film fluid dynamics near Marx's line.

Extensive studies have explored the dynamics of the ocular surface fluid, though theoretical investigations are typically limited to the use of the lubrication approximation, which is not guaranteed to be uniformly valid a-priori throughout the tear meniscus. However, resolving tear film behaviour within the meniscus and especially its apices is required to characterise the flow dynamics where the tear film is especially thin, and thus most susceptible to evaporatively induced hyperosmolarity and subsequent epithelial damage. Hence, we have explored the accuracy of the standard lubrication approximation for the tear film by explicit comparisons with the 2D Navier-Stokes model, considering both stationary and moving eyelids. Our results demonstrate that the lubrication model is qualitatively accurate except in the vicinity of the eyelids. In particular, and in contrast to lubrication theory, the solution of the full Navier-Stokes equations predict a distinct absence of fluid flow, and thus convective mixing in the region adjacent to the tear film contact line. These observations not only support emergent hypotheses concerning the formation of Marx's line, a region of epithelial cell staining adjacent to the contact line on the eyelid, but also enhance our understanding of the pathophysiological consequences of the flow profile near the tear film contact line.

Predicting the safety and efficacy of buffer therapy to raise tumour pHe: an integrative modelling study.

BACKGROUND: Clinical positron emission tomography imaging has demonstrated the vast majority of human cancers exhibit significantly increased glucose metabolism when compared with adjacent normal tissue, resulting in an acidic tumour microenvironment. Recent studies demonstrated reducing this acidity through systemic buffers significantly inhibits development and growth of metastases in mouse xenografts. METHODS: We apply and extend a previously developed mathematical model of blood and tumour buffering to examine the impact of oral administration of bicarbonate buffer in mice, and the potential impact in humans. We recapitulate the experimentally observed tumour pHe effect of buffer therapy, testing a model prediction in vivo in mice. We parameterise the model to humans to determine the translational safety and efficacy, and predict patient subgroups who could have enhanced treatment response, and the most promising combination or alternative buffer therapies. RESULTS: The model predicts a previously unseen potentially dangerous elevation in blood pHe resulting from bicarbonate therapy in mice, which is confirmed by our in vivo experiments. Simulations predict limited efficacy of bicarbonate, especially in humans with more aggressive cancers. We predict buffer therapy would be most effectual: in elderly patients or individuals with renal impairments; in combination with proton production inhibitors (such as dichloroacetate), renal glomular filtration rate inhibitors (such as non-steroidal anti-inflammatory drugs and angiotensin-converting enzyme inhibitors), or with an alternative buffer reagent possessing an optimal pK of 7.1-7.2. CONCLUSION: Our mathematical model confirms bicarbonate acts as an effective agent to raise tumour pHe, but potentially induces metabolic alkalosis at the high doses necessary for tumour pHe normalisation. We predict use in elderly patients or in combination with proton production inhibitors or buffers with a pK of 7.1-7.2 is most promising.

Modelling a tethered mammalian sperm cell undergoing hyperactivation.

The beat patterns of mammalian sperm flagella can be categorised into two different types. The first involves symmetric waves propagating down the flagellum with a net linear propulsion of the sperm cell. The second, hyperactive, waveform is classified by vigorous asymmetric waves of higher amplitude, lower wavenumber and frequency propagating down the flagellum resulting in highly curved trajectories. The latter beat pattern is part of the capacitation process whereby sperm prepare for the prospective penetration of the zona pellucida and fusion with the egg. Hyperactivation is often observed to initiate as sperm escape from epithelial and ciliary bindings formed within the isthmic regions of the female oviducts, leading to a conjecture in the literature that this waveform is mechanically important for sperm escape. Hence, we explore the mechanical effects of hyperactivation on a tethered sperm, focussing on a Newtonian fluid. Using a resistive force theory model we demonstrate that hyperactivation can indeed generate forces that pull the sperm away from a tethering point and consequently a hyperactivated sperm cell bound to an epithelial surface need not always be pushed by its flagellum. More generally, directions of the forces generated by tethered flagella are insensitive to reductions in beat frequency and the detailed flagellar responses depend on the nature of the binding at the tethering point. Furthermore, waveform asymmetry and amplitude increases enhance the tendency for a tethered flagellum to start tugging on its binding. The same is generally predicted to be true for reductions in the wavenumber of the flagellum beat, but not universally so, emphasising the dynamical complexity of flagellar force generation. Finally, qualitative observations drawn from experimental data of human sperm bound to excised female reproductive tract are also presented and are found to be consistent with the theoretical predictions.

Re-evaluating the use of Voronoi Tessellations in the assessment of oxygen supply from capillaries in muscle.

The ability to characterise capillary supply plays a key role in developing effective therapeutic interventions for numerous pathological conditions, such as capillary loss in skeletal or cardiac muscle. However, quantifying capillary supply is fraught with difficulties. Averaged measures such as capillary density or mean inter-capillary distance cannot account for the local geometry of the underlying capillary distribution, and thus can only highlight a tissue wide, global hypoxia. Detailed tissue geometry, such as muscle fibre size, has been incorporated into indices of capillary supply by considering the distribution of Voronoi tessellations generated from capillary locations in a plane perpendicular to muscle fibre orientation, implicitly assuming that each Voronoi polygon represents the area of supply of its enclosed capillary. Using a modelling framework to assess the capillary supply capacity under maximal sustainable conditions in muscle, we theoretically demonstrate that Voronoi tessellations often provide an accurate representation of the regions supplied by each capillary. However, we highlight that this use of Voronoi tessellations is inappropriate and inaccurate in the presence of extensive capillary rarefaction and pathological variations in oxygen tension of different capillaries. In such cases, oxygen flux trapping regions are developed to provide a more general representation of the capillary supply regions, in particular incorporating the additional influences of heterogeneity that are absent in the consideration of Voronoi tessellations.

Coupling fluid and solute dynamics within the ocular surface tear film: a modelling study of black line osmolarity.

We present a mathematical model describing the spatial distribution of tear film osmolarity across the ocular surface of a human eye during one blink cycle, incorporating detailed fluid and solute dynamics. Based on the lubrication approximation, our model comprises three coupled equations tracking the depth of the aqueous layer of the tear film, the concentration of the polar lipid, and the concentration of physiological salts contained in the aqueous layer. Diffusive boundary layers in the salt concentration occur at the thinnest regions of the tear film, the black lines. Thus, despite large Peclet numbers, diffusion ameliorates osmolarity around the black lines, but nonetheless is insufficient to eliminate the build-up of solute in these regions. More generally, a heterogeneous distribution of solute concentration is predicted across the ocular surface, indicating that measurements of lower meniscus osmolarity are not globally representative, especially in the presence of dry eye. Vertical saccadic eyelid motion can reduce osmolarity at the lower black line, raising the prospect that select eyeball motions more generally can assist in alleviating tear film hyperosmolarity. Finally, our results indicate that measured evaporative rates will induce excessive hyperosmolarity at the black lines, even for the healthy eye. This suggests that further evaporative retardation at the black lines, for instance due to the cellular glycocalyx at the ocular surface or increasing concentrations of mucus, will be important for controlling hyperosmolarity as the black line thins.

Comment on the article by J. Elgeti, U. B. Kaupp, and G. Gompper: hydrodynamics of sperm cells near surfaces.

A recent study by Elgeti et al. used multiparticle collision dynamics to simulate a long-standing problem: the approach of sperm to surfaces, and subsequent accumulation. The authors highlight differences in their predictions with those of the earlier Stokes flow simulations of Smith et al. attributing the differences to methodological flaws in the earlier article. In this Comment, we discuss the criticisms leveled in detail, and review some recently published work that shows how species-specific details of cell morphology provides a more likely explanation for the differing predictions of the two studies. We also highlight experimental work that supports the study of Smith et al.

Tumour angiogenesis: the gap between theory and experiments.

A common experimental technique for viewing in vivo angiogenesis utilises tumours implanted into a test animal cornea. The cornea is avascular but the tumour promotes vascularisation from the limbus and the new blood vessels can be readily observed through the transparent cornea. Many of the early mathematical models for tumour angiogenesis used this scenario as their experimental template and as such assumed that there is a large gap, of the order of 2mm, between the tumour and neighbouring vasculature at the onset of angiogenesis. In this work we consider whether the assumption that there is a significant gap between the tumour and neighbouring vasculature is unique to intra-cornea tumour implants, or whether this characterises avascular tumour growth more generally. To do this we utilise a simple scaling argument, derive a multi-compartment model for tumour growth, and consider in vivo images. This analysis demonstrates that the corneal implant experiments and the corresponding mathematical models cannot generally be applied to a clinical setting.

The dynamics of Turing patterns for morphogen-regulated growing domains with cellular response delays.

Since its conception in 1952, the Turing paradigm for pattern formation has been the subject of numerous theoretical investigations. Experimentally, this mechanism was first demonstrated in chemical reactions over 20 years ago and, more recently, several examples of biological self-organisation have also been implicated as Turing systems. One criticism of the Turing model is its lack of robustness, not only with respect to fluctuations in the initial conditions, but also with respect to the inclusion of delays in critical feedback processes such as gene expression. In this work we investigate the possibilities for Turing patterns on growing domains where the morphogens additionally regulate domain growth, incorporating delays in the feedback between signalling and domain growth, as well as gene expression. We present results for the proto-typical Schnakenberg and Gierer-Meinhardt systems: exploring the dynamics of these systems suggests a reconsideration of the basic Turing mechanism for pattern formation on morphogen-regulated growing domains as well as highlighting when feedback delays on domain growth are important for pattern formation.

Aberrant behaviours of reaction diffusion self-organisation models on growing domains in the presence of gene expression time delays.

Turing's pattern formation mechanism exhibits sensitivity to the details of the initial conditions suggesting that, in isolation, it cannot robustly generate pattern within noisy biological environments. Nonetheless, secondary aspects of developmental self-organisation, such as a growing domain, have been shown to ameliorate this aberrant model behaviour. Furthermore, while in-situ hybridisation reveals the presence of gene expression in developmental processes, the influence of such dynamics on Turing's model has received limited attention. Here, we novelly focus on the Gierer-Meinhardt reaction diffusion system considering delays due the time taken for gene expression, while incorporating a number of different domain growth profiles to further explore the influence and interplay of domain growth and gene expression on Turing's mechanism. We find extensive pathological model behaviour, exhibiting one or more of the following: temporal oscillations with no spatial structure, a failure of the Turing instability and an extreme sensitivity to the initial conditions, the growth profile and the duration of gene expression. This deviant behaviour is even more severe than observed in previous studies of Schnakenberg kinetics on exponentially growing domains in the presence of gene expression (Gaffney and Monk in Bull. Math. Biol. 68:99-130, 2006). Our results emphasise that gene expression dynamics induce unrealistic behaviour in Turing's model for multiple choices of kinetics and thus such aberrant modelling predictions are likely to be generic. They also highlight that domain growth can no longer ameliorate the excessive sensitivity of Turing's mechanism in the presence of gene expression time delays. The above, extensive, pathologies suggest that, in the presence of gene expression, Turing's mechanism would generally require a novel and extensive secondary mechanism to control reaction diffusion patterning.

The influence of gene expression time delays on Gierer-Meinhardt pattern formation systems.

There are numerous examples of morphogen gradients controlling long range signalling in developmental and cellular systems. The prospect of two such interacting morphogens instigating long range self-organisation in biological systems via a Turing bifurcation has been explored, postulated, or implicated in the context of numerous developmental processes. However, modelling investigations of cellular systems typically neglect the influence of gene expression on such dynamics, even though transcription and translation are observed to be important in morphogenetic systems. In particular, the influence of gene expression on a large class of Turing bifurcation models, namely those with pure kinetics such as the Gierer-Meinhardt system, is unexplored. Our investigations demonstrate that the behaviour of the Gierer-Meinhardt model profoundly changes on the inclusion of gene expression dynamics and is sensitive to the sub-cellular details of gene expression. Features such as concentration blow up, morphogen oscillations and radical sensitivities to the duration of gene expression are observed and, at best, severely restrict the possible parameter spaces for feasible biological behaviour. These results also indicate that the behaviour of Turing pattern formation systems on the inclusion of gene expression time delays may provide a means of distinguishing between possible forms of interaction kinetics. Finally, this study also emphasises that sub-cellular and gene expression dynamics should not be simply neglected in models of long range biological pattern formation via morphogens.

On the modelling of biological patterns with mechanochemical models: Insights from analysis and computation.

The diversity of biological form is generated by a relatively small number of underlying mechanisms. Consequently, mathematical and computational modelling can, and does, provide insight into how cellular level interactions ultimately give rise to higher level structure. Given cells respond to mechanical stimuli, it is therefore important to consider the effects of these responses within biological self-organisation models. Here, we consider the self-organisation properties of a mechanochemical model previously developed by three of the authors in Acta Biomater. 4, 613-621 (2008), which is capable of reproducing the behaviour of a population of cells cultured on an elastic substrate in response to a variety of stimuli. In particular, we examine the conditions under which stable spatial patterns can emerge with this model, focusing on the influence of mechanical stimuli and the interplay of non-local phenomena. To this end, we have performed a linear stability analysis and numerical simulations based on a mixed finite element formulation, which have allowed us to study the dynamical behaviour of the system in terms of the qualitative shape of the dispersion relation. We show that the consideration of mechanotaxis, namely changes in migration speeds and directions in response to mechanical stimuli alters the conditions for pattern formation in a singular manner. Furthermore without non-local effects, responses to mechanical stimuli are observed to result in dispersion relations with positive growth rates at arbitrarily large wavenumbers, in turn yielding heterogeneity at the cellular level in model predictions. This highlights the sensitivity and necessity of non-local effects in mechanically influenced biological pattern formation models and the ultimate failure of the continuum approximation in their absence.

Bend propagation in the flagella of migrating human sperm, and its modulation by viscosity.

A pre-requisite for sexual reproduction is successful unification of the male and female gametes; in externally-fertilising echinoderms the male gamete is brought into close proximity to the female gamete through chemotaxis, the associated signalling and flagellar beat changes being elegantly characterised in several species. In the human, sperm traverse a relatively high-viscosity mucus coating the tract surfaces, there being a tantalising possible role for chemotaxis. To understand human sperm migration and guidance, studies must therefore employ similar viscous in vitro environments. High frame rate digital imaging is used for the first time to characterise the flagellar movement of migrating sperm in low and high viscosities. While qualitative features have been reported previously, we show in precise spatial and temporal detail waveform evolution along the flagellum. In low viscosity the flagellum continuously moves out of the focal plane, compromising the measurement of true curvature, nonetheless the presence of torsion can be inferred. In high viscosities curvature can be accurately determined and we show how waves propagate at approximately constant speed. Progressing waves increase in curvature approximately linearly except for a sharper increase over a distance approximately 20-27 microm from the head/midpiece junction. Curvature modulation, likely influenced by the outer dense fibres, creates the characteristic waveforms of high viscosity swimming, with remarkably effective cell progression against greatly increased resistance, even in high viscosity liquids. Assessment of motility in physiological viscosities will be essential in future basic research, studies of chemotaxis and novel diagnostics.

Modelling mucociliary clearance.

Mathematical modelling of the fluid mechanics of mucociliary clearance (MCC) is reviewed and future challenges for researchers are discussed. The morphology of the bronchial and tracheal airway surface liquid (ASL) and ciliated epithelium are briefly introduced. The cilia beat cycle, beat frequency and metachronal coordination are described, along with the rheology of the mucous layer. Theoretical modelling of MCC from the late 1960s onwards is reviewed, and distinctions between 'phenomenological', 'slender body theory' and recent 'fluid-structure interaction' models are explained. The ASL consists of two layers, an overlying mucous layer and underlying watery periciliary layer (PCL) which bathes the cilia. Previous models have predicted very little transport of fluid in the PCL compared with the mucous layer. Fluorescent tracer transport experiments on human airway cultures conducted by Matsui et al. [Matsui, H., Randell, S.H., Peretti, S.W., Davis, C.W., Boucher, R.C., 1998. Coordinated clearance of periciliary liquid and mucus from airway surfaces. J. Clin. Invest. 102 (6), 1125-1131] apparently showed equal transport in both the PCL and mucous layer. Recent attempts to resolve this discrepancy by the present authors are reviewed, along with associated modelling findings. These findings have suggested new insights into the interaction of cilia with mucus due to pressure gradients associated with the flat PCL/mucus interface. This phenomenon complements previously known mechanisms for ciliary propulsion. Modelling results are related to clinical findings, in particular the increased MCC observed in patients with pseudohypoaldosteronism. Recent important advances by several groups in modelling the fluid-structure interaction by which the cilia movement and fluid transport emerge from specification of internal mechanics, viscous and elastic forces are reviewed. Finally, we discuss the limitations of existing work, and the challenges for the next generation of models, which may provide further insight into this complex and vital system.

A mass action model of a Fibroblast Growth Factor signaling pathway and its simplification.

We consider a kinetic law of mass action model for Fibroblast Growth Factor (FGF) signaling, focusing on the induction of the RAS-MAP kinase pathway via GRB2 binding. Our biologically simple model suffers a combinatorial explosion in the number of differential equations required to simulate the system. In addition to numerically solving the full model, we show that it can be accurately simplified. This requires combining matched asymptotics, the quasi-steady state hypothesis, and the fact subsets of the equations decouple asymptotically. Both the full and simplified models reproduce the qualitative dynamics observed experimentally and in previous stochastic models. The simplified model also elucidates both the qualitative features of GRB2 binding and the complex relationship between SHP2 levels, the rate SHP2 induces dephosphorylation and levels of bound GRB2. In addition to providing insight into the important and redundant features of FGF signaling, such work further highlights the usefulness of numerous simplification techniques in the study of mass action models of signal transduction, as also illustrated recently by Borisov and co-workers (Borisov et al. in Biophys. J. 89, 951-966, 2005, Biosystems 83, 152-166, 2006; Kiyatkin et al. in J. Biol. Chem. 281, 19925-19938, 2006). These developments will facilitate the construction of tractable models of FGF signaling, incorporating further biological realism, such as spatial effects or realistic binding stoichiometries, despite a more severe combinatorial explosion associated with the latter.