PDE/ODE modeling and simulation to determine the role of diffusion in longterm and range cellular signaling
 Elfriede Friedmann^{1}Email author
https://doi.org/10.1186/s1362801500248
© Friedmann. 2015
Received: 5 November 2014
Accepted: 21 September 2015
Published: 14 October 2015
Abstract
Background
We study the relevance of diffusion for the dynamics of signaling pathways. Mathematical modeling of cellular diffusion leads to a coupled system of differential equations with Robin boundary conditions which requires a substantial knowledge in mathematical theory. Using our new developed analytical and numerical techniques together with modern experiments, we analyze and quantify various types of diffusive effects in intra and intercellular signaling. The complexity of these models necessitates suitable numerical methods to perform the simulations precisely and within an acceptable period of time.
Methods
The numerical methods comprise a Galerkin finite element space discretization, an adaptive time stepping scheme and either an iterative operator splitting method or fully coupled multilevel algorithm as solver.
Results
The simulation outcome allows us to analyze different biological aspects. On the scale of a single cell, we showed the high cytoplasmic concentration gradients in irregular geometries. We found an 11 % maximum relative total STAT5concentration variation in a fibroblast and a 70 % maximum relative pSTAT5concentration variation in a fibroblast with an irregular cell shape. For pSMAD2 the maximum relative variation was 18 % in a hepatocyte with a box shape and 70 % in an irregular geometry. This result can be also obtained in a cell with a box shape if the molecules diffuse slowly (with D=1 μm^{2}/s instead of D=15 μm^{2}/s). On a scale of cell system in the lymph node, our simulations showed an inhomogeneous IL2 pattern with an amount over three orders of magnitude (10^{−3}−1 pM) and high gradients in face of its fast diffusivity. We observed that 20 out of 125 cells were activated after 9 h and 33 in the steady state. Our insilico experiments showed that the insertion of 31 regulatory T cells in our cell system can completely downregulate the signal.
Conclusions
We quantify the concentration gradients evolving from the diffusion of the molecules in several signaling pathways. For intracellular signaling pathways with nuclear accumulation the size of cytoplasmic gradients does not indicate the change in gene expression which has to be analyzed separately in future. For intercellular signaling the high cytokine concentration gradients play an essential role in the regulation of the molecular mechanism of the immune response. Furthermore, our simulation results can give the information on which signaling pathway diffusion may play a role. We conclude that a PDE model has to be considered for cells with an irregular shape or for slow diffusing molecules. Also the high gradients inside a cell or in a cell system can play an essential role in the regulation of the molecular mechanisms.
Keywords
Systems of coupled differential equations Reactiondiffusion systems Numerical simulation Cellular signalingBackground
A single cell, the smallest unit of life, alters, learns, adapts, specializes and differentiates. One of the causes of the diversification of a cell is the concentration distribution of signaling molecules which trigger signaling pathways. Changes in signaling pathways and in the amount of molecule concentration induce different reactions in the specific genes. These genes are controlled by proteins (transcription factors) and play an important role on the differentiation. They, in turn, control the concentration as well as the moment and duration of activation. To understand these biological processes it requires a detailed analysis of the signal transduction processes and their interplay. These complex processes are commonly described by mathematical models to generate experimentally testable hypotheses. With the available data different mathematical models can be derived. Initially in signal transduction, simple models consisting mainly of cascades of events were developed. Then, more general models followed to describe mechanisms such as pattern formation and regulation of immune response. A review on the state of the mathematical models in signal transduction can be found in [1]. In most cases, ordinary differential equations (ODE) are used to identify the role of specific components of the pathway and to determine the dynamics and outcome predictions. Partial differential equations (PDE) were used to analyze pattern formation for example stripe formation in zebras or fish [2], and cell structure properties like deformability, cell polarity [3] or cell migration [4].
In this paper, we use the reaction diffusion system for cellular signaling pathways to find the spatial distribution of the signaling molecules which causes the diversification of cells and their actions. Our contribution is to include the spatial aspect (size and geometry of the cell) of the databased signal transduction models considered by our collaborators. Also we develop the mathematical models that describe the key dynamic properties and predict strategies for intervention. We extend a simple ODEsystem to a more complicated PDEsystem of reactiondiffusion equations. In PDE we take into account the traveling of the signal molecules from the outer membrane through the cytoplasm, thus forming a concentration gradient within the cell. The focus remains on the spatial distributions of the important signaling molecules whereas other processes can be replaced by specific measurements (e.g. the phosphorylated JAK (pJAK) concentration can be determined and replace the whole receptor model) or by ODE if we assume that the local domain (nucleus) of the interplay of the concentrations is considered to be well mixed (i.e. all processes take place everywhere).
The modeling leads to a coupled PDE/ODE system with Robin type boundary conditions. This nonstandard model demands a substantial knowledge in mathematical theory. It needs the additional development of both analytical tools and numerical algorithms for its analysis, simulation and evaluation, which explains why the model is not commonly considered in other fields. Moreover, additional changes can occur on the dynamics of the pathway when a diffusion is included as shown in [5] for the CalvinCycle. The diffusion of the molecules can cause a gradient in their concentration which often plays an important role such as the one of RanGTPImportin β. This gradient affects the stabilization of microtubules and the formation of the mitotic spindle [6, 7].
We are interested in biological processes regulated by diffusion of molecules over great distances (many cell distances) and a long period of time (over 200 min in the intracellular and 30 h in the intercellular pathways). We start with the classical diffusion model based on Fick’s law [8]. Given the diffusion coefficient the velocity of the diffusion depends on the nature of the medium, the size of the particles and temperature. In order to describe the complicated processes within a cell we consider a normal diffusion, an aided diffusion (such as energycontrolled movement or active transport with help of transport molecules) and a retarded diffusion (due to the binding process on other molecules). This leads to a nonuniform spatial distribution of molecules. Unfortunately, the current microscopy methods are not fully capable of measuring all dynamic parameters. Thus, diffusion maps inside the cell based on the experimental data are not yet available. In this paper we use the diffusion coefficient measured by our collaborators via Fluorescence Recovery After Photobleaching (FRAP) and Fluorescence Correlation Spectroscopy (FCS). In our models the diffusion coefficient has hence an extended physical meaning (it is a macroscopic approximation of the cytoplasmic processes) and can deviate from the one used in pure biomolecular simulations.
Our databased modeling simplifies the system of equations because the certain biological phenomena can be replaced by a measured quantity (e.g. the receptor model in the JAK2/STAT5 pathway was skipped due to the measurements of the concentration of the phosphorylated JAK2molecules) [9]. In the models we use the highquality quantitative data from our collaborators [10]. This simplification is necessary for the theoretical treatment of the model and the numerics in three dimensions. The derived models have a structure allowing the mathematical theory to be suitable for several applications. For example, there is a particle exchange between two domains, one wellmixed domain (ODE model) and an adjacent nonwell mixed domain (PDE model). Reactions occur in the wellmixed domain and reaction and diffusion occur in the nonwell mixed domain.
In this paper we discuss whether the diffusion of signaling molecules plays a role in the signaling of intra and intercellular systems. To illustrate this aspect we consider the models of the wellknown pathways such as JAK2/STAT5, SMAD and IL2. These nonstandard models need the additional development on both analytical tools and numerical algorithms for further analysis, simulation and evaluation. The new analytical and numerical tools are presented in [11, 12]. The simulation results of the JAK/STAT pathway are presented for the first time in this work. We quantify the influence of the cell geometry and diffusivity on the dynamics of the signal. We also quantify the concentration gradients during the interplay between diffusion and degradation in order to discuss its effect to the cell fate. The biological interpretation of the simulation results are presented in a compact form to provide an overview in which pathways and cells the diffusion may play a role.
The paper has five main sections. After an introduction to the subject in the first section, the section Methods presents two of our mathematical tools: modeling and numerical simulations. We present the general mathematical models which differ from the commonly used models for these applications. We will give a short overview of their mathematical analysis and numerical calculations to underline the importance of the theoretical background. We present our numerical methods based on finite elements which are used to simulate the modeled equations. We use a different method to solve the stronger nonlinear coupling in the intercellular model. Therefore, each of the section Results and Discussion comprises two subsections where we examine models with intra and intercellular diffusion separately. An analysis of the numerical methods used will be subject of a forthcoming paper. The section Discussion contains the observations and the biological interpretation of the numerical results.
Methods
Model setup
Since the process inside a cell is highly complex all the existing models are only an attempt of an approximation. Membrane receptors receive extracellular signals from each cell in form of protein concentrations. The signals are processed, encoded and transferred to the nucleus where a further development is decided. This signal processing from the cell membrane to the nucleus occurs via activation (phosphorylation) and spatial relocation of the components of the signaling pathways. The cytoplasmic proteins are recruited to the cell membrane where the activation takes place. The activated molecules dissociate from the receptor, dimerize, trimerize or form other complexes and move to the nucleus to regulate the transcription of target genes. Then, the molecules are deactivated and are shuttled back to the cytoplasm where the whole process restarts until the cell cycle in the considered pathway is regulated. The three factors that influence the process are the fast diffusion rate, the activation rate and the boundary conditions at the interface between cytoplasm and nucleus which determine a nuclear accumulation of the activated molecules. A steady state is reached by the continuous molecule shuttling.
In [6] it shows that the spatial separation of activation and deactivation mechanisms can result in steep phosphoprotein gradients. For spherical cells, a single concentration and cytosolic phosphatase, diffusion must be taken into account if D/L ^{2}<<k, where L describes the traveled distance and k the dissociation rate of the molecules. In [12] it shows that this estimation can be used for our model in the case of cells with a regular shape. For general cases and nonlinear models this estimation can not be applied, and we need three dimensional simulations for each considered pathway. The cell fate is not only determined by the cell itself but also by surrounding cells (celltocell communication) and molecules and a complex interplay between them. We will present the models for signaling pathways inside a cell (intracellular signaling pathways) and between cells (intercellular signaling pathways) that appear to have the same mathematical structure.
The coupled PDE/ODE reaction diffusion model
for 1≤p<∞, and \(\,\u\_{L^{\infty }(0,T;X)} := \text {ess sup}_{(0,T)} \,\u(t)\<\infty \). We denote by (·,·)_{ Ω } and (·,·)_{ Γ } the usual L ^{2} scalar products over the domain Ω and a part Γ⊂∂ Ω of its boundary, respectively, and by ∥·∥_{ Ω } and ∥·∥_{ Γ } the corresponding norms.
The detailed model description for the JAK2/STAT5 signaling pathway can be found in [9], the SMAD pathway in [12] and the IL2 signaling in [13–16].
The mathematical analysis and numerical simulation of such systems can not be performed with a standard theory. The methods depend very much on the variety factors such as structure of the equations, nonlinearities, coupling, domain and more. For a linear model for JAK2/STAT5, we showed the wellposedness and the properties of the solution for reliable numerical simulation [11]. In [17] analytical investigations are performed for a simplified situation of two cells in two dimensions. In [14, 17] numerical simulations were performed in two dimension by using Finite Differences (FD). Nevertheless, the two dimensional results represent an artificial reality, e.g. the 2D cells are infinitelylong cylinders instead of spheres. The cytokine Interleukin 2 (IL2) can spread only in two directions causing the change in its amount and gradient. For a better understanding of the process, realistic simulations in three dimension are necessary [13].
In the next section we present the numerical simulations of (2)–(6).
Numerical simulations
General description of the methods

Galerkin space discretization by Finite Elements (FE)

segregating solution approach [12] or a fully coupled multilevel algorithm as solver [13]

multilevel preconditioner with a BlockGaussSeidel scheme as smoother consisting of damped Jacobi iterations on the PDE part (presented in a forthcoming paper)

error control techniques with the Dual Weighted Residual (DWR) method for space and time adaptive discretization [20, 21].
A segregating solution approach is a splitting solution approach. It is often used when the restriction on accuracy can be relaxed in order to allow an easier numerical treatment of complicated problems. Such an approach makes it possible to reuse the existing solvers and is widely used in numerical methods for coupled systems. The equations are discretized via the Rothe Method, first on temporal discretization by applying backward Euler or a similar solver of second order, then on spatial discretization via Finite Elements. The two parts are solved sequentially with splitted operators which keep the storage space low. In case of a strong coupling, it requires a very small time step which leads to a long computing time. Therefore, we developed an adaptive and fully coupled solver for the intercellular model.
Because the main focus of this paper is devoted to biological diffusion, we will keep the mathematical aspects short. Interested readers are recommended to look trough the literature cited.
Results
Intracellular diffusion
The model for the different cell types requires different input parameters. Some parameters were measured by specific experimental methods and others were estimated [9]. Fibroblasts can have various shapes often depending on the current function or tissue concentration [12]. This has to be considered in modeling as well as in the experiments when a certain parameter is measured.
One of the advantages of numerical simulations (in silico experiments) is that we gain additional information about the processes modeled. The outcome of the simulation can be adapted to meet the requirement of the experiment. For this pathway both states of the molecules (non and activated state, u _{0}+u _{1}) must be marked and measured together. The quantity of each specie or other quantities of interest that can not be evaluated in the experiment can be determined from the simulation.
Result 1.
For a perfectly spherical cell such as the CFUE cell with a large nucleus compared to the amount of cytoplasm, we observed a homogeneous distribution of the signaling molecules due to the fast diffusion (Fig. 2). For this geometry we observed no deviation between the ODE and the PDE model.
Result 2.
For the modeled fibroblast we observed a maximum relative variation of 11 % of the total STATconcentration. The concentration in the cytoplasm was about 1.5 molecules/μm^{3} higher compared to the outcome of the pure ODE model due to the diffusion of the molecules into the nucleus. The activation process takes a place at the cell receptors located at the outer membrane which concludes that the main contribution to this concentration variation comes from the activated STAT5concentration. Comparing the results of the ODE model with the PDE model the largest difference in the activated STAT5concentration was 0.7 molecules/μm^{3} for the modeled fibroblast (Fig. 3). This amount occurred 10 min after starting the activation process. During 200 min of simulation a steady state was observed in the last 80 min.
Result 3.
The simulation result for the reconstructed geometry lies in between the results of the CFUE cell and the modeled fibroblast. Here, we obtained a maximum relative variation of 10 % of the activated STATconcentration, a variation of 2−3.5 % of the nonactivated STATconcentration and a concentration difference of 0.07 molecules/μm^{3} between the ODE and PDE model (Fig. 4).
Result 4.
The greatest relative concentration variation was observed in the artificial cell assembled with the two extensions: 40 % for the unphosphorylated STAT5molecules and 70 % for the phosphorylated (Fig. 5). The deviation in the pSTAT5concentration between the two models (ODEPDE) was five times greater (0.35 molecules/μm^{3}) than in the same cell without extensions.
Remark.
We choose the coefficients \(\tilde {D}_{\textit {xx}}, \tilde {D}_{\textit {yy}}\) and \(\tilde {D}_{\textit {zz}}\) so that we have higher diffusion towards the nucleus \(\tilde {D}_{\textit {yy}}>> D\) and smaller diffusion elsewhere \(\tilde {D}_{\textit {xx}}= \tilde {D}_{\textit {zz}}<<D\). Also we want the same trace in the diffusion tensors, i.e. \(\tilde {D}_{\textit {xx}}+\tilde {D}_{\textit {yy}}+\tilde {D}_{\textit {zz}}=3D\). The results are presented in Fig. 6:
Result 5.
For spatial inhomogeneous diffusion of the molecules we obtained the same maximum relative variation of the total STATconcentration (1.5 molecules/μm^{3}). However, the deviation between the ODE and PDE model was six times greater after 10 min of activation (0.42 molecules/μm^{3}) and three times greater in the steady state (0.23 molecules/μm^{3}) than in the cell with isotropic diffusion. This coincides with the fact that some molecules reach the nucleus faster but others have a longer sojourn time in the cytoplasm.
The hepatocyte model represents a cell from the liver. It has a regular boxshape due to the high cell density in the tissue. The reconstructed cell is flatter and elongated due to microscopy reconstruction from a cell culture with a lessdense cell concentration [12]. In the artificial cell we added two cytoplasmic extensions to the modeled geometry to see any possible effect of these cytoplasmic structures which commonly appear in lessdense cell cultures.
Result 6.
We observed a much greater maximum relative concentration variation in the activated SMAD concentration than in the activated STAT concentration in the JAK2/STAT5 pathway. For the activated SMAD2 we obtained 18 % in the box shaped cell, 38 % in the reconstructed and 50 % in the artificial cell [12].
Result 7.
A visible difference in the nuclear SMADtrimer concentration between the ODE and PDE model was seen only in the artificial cell with extensions.
Result 8.
For the modeled hepatocyte with a regular form we observed no visible difference in the nuclear SMADtrimer concentration between the ODE and PDE model except for slower diffusing molecules.
Different diffusion coefficients give rise to different concentration distributions:
Result 9.
We will give a more detailed discussion about the size of cytosolic gradients and their effect to the development of the cell in the next section, Discussion about biological diffusion.
Intercellular diffusion
Our model captures the first 30 h of IL2 signaling, the initial phase after antigen stimulation where the cells are primed for proliferation but have not yet entered initiated cell division. Biological details can be found in [14–16].
Result 10.
Result 11.
Result 12.
The Fig. 14 shows that only 37 cells were activated in immediate vicinity of the secreting cell (marked with yellow) by a secretion rate of 10^{6} molecules/h which confirms the strong localized uptake of IL2. Our calculations showed that with a higher secretion rate (1.5×10^{6} molecules/h) 123 cells were activated.
In the next section we will discuss the importance of the modeled biological diffusion for the intra and as well for the intercellular signaling pathways.
Discussion about biological diffusion
Intracellular diffusion
Our simulations give us the potential to analyze the effect of diffusion of signaling molecules. The signal to the nucleus may be delayed or modulated by diffusive processes which causes changes in gene expression. The possibility of an impact of the changed signal on the further development of the cell has to be analyzed separately in each model.
In case of the SMAD signaling pathway we investigate further on the role of diffusion. In addition to the nuclear accumulation, our intracellular models exhibit mass conservation. The molecules change the compartments (cytoplasm and nucleus), their state (nonactivated to activated and vice versa) and interact with each other while their number is preserved. The deviation of the concentrations (pSMAD2, trimer) in the cytoplasm is balanced by a deviation of the corresponding concentrations in the nucleus. Due to the nuclear accumulation any remarkable cytoplasmic concentration deviation is slightly small in the nucleus for cells with regular shape. In this case enough transcription factors may be available in the nucleus and do not have to compete for the binding sides of the DNA except for the case where the molecules diffuses slower than D=1 μ m ^{2}/s. For a diffusion coefficient D=15 μ m ^{2}/s we observed in Fig. 9 no change in the nuclear concentration in the PDE model compared to the ODE model and for D=1 μ m ^{2}/s a small change with a possible effect to the signal which we will discuss later.
Result 13.
Our simulation results showed that in signaling pathways with nuclear accumulation the size of cytoplasmic gradients does not indicate the change in gene expression. For cells with regular shape the effect of cell geometry to the signal can be neglected. The concentration difference in the nucleus where gene expression occurs is small and the effect to the further development of the cell is not known yet.
Result 14.
In a cell with a regular shape the slow diffusing molecules (D<1 μ m ^{2}/s) cause the change in the nuclear concentration.
Result 15.
For cells with irregular shape the effect of diffusion is greater. The cell geometry has greater influence on the change in the nuclear concentration than the diffusion velocity.
Intercellular diffusion
In contrast to the intracellular signaling pathways discussed above, we will see that for the intercellular IL2 signaling pathway high gradients are crucial for the effective signal.
Upon stimulation by antigen in the lymph node, the proliferation and the differentiation of T lymphocytes are tightly regulated by the interplay of regulatory T cells and T helper cells. These T helper cells recognize an antigen of an antigen presenting cell and then activate responding T cells as a part of the immune response. Regulatory T cells instead suppress the activation of the responding T cells. Understanding the role of regulatory T cells is important for the treatment of autoimmune diseases and cancer. The success of organ transplantation and cancer immunotherapy are directly linked to the suppressing activity of regulatory T cells. The experimental, theoretical [14] and numerical studies [17] have identified secretion and uptake of IL2 as a possible mechanism mediating immune suppression by regulatory T cells.
In [25] it was found that the cytokines are secreted in a polarized way at the immunological synapse preferred in a very narrow space nearby the contact with the antigen presenting cell. In particular, it is not understood within which spatial range cytokines can signal. It was speculated that the range of influence of the cytokine remains restricted to these two contacting cells. On the other side, it was shown that all cells in the lymph node can sense the cytokine IL4 [26]. This gives the evidence to a more global distribution of cytokines.
Cytokine concentrations measured by ELISA studies in T cell cultures are typically in the pM range [27]. In this regime, cytokine molecules can reach their targets only after the diffusion over large distance as seen from the molecular scale. We assume that the time for IL2 molecules to diffuse towards a T cell is short. In contrast, the average time to reach a receptor at the cell surface is long because naive T cells express only small number of cytokine receptors, which is not detectable in experiments. Thus, IL2 concentrations measured by ELISA are too low to generate reliable signals. It has been reported that IL2 is subject to huge gradients with much higher concentration at the surface of T cells [28]. To investigate such inhomogeneities, it is necessary to consider the spatiotemporal dynamics of cytokine signals. More details are given in [13].
We determine the number of activated cells (Fig. 12), the amount of IL2 in the extracellular space (Fig. 13), the range of action of IL2 (Fig. 14) and the required number of regulatory cells to downregulate the signal (Fig. 12). In our calculation we observe a variety of IL2 amount over three orders of magnitude over the entire simulation time (Fig. 13 left). Thus, for activation the cells require both a transient strong and a stable weak IL2 signal. The competition for IL2 in the lymph node among regulatory T cells, responding T cells and T helper cells is thus a very local process. The position of the secretory cells decides which cell will be activated in the absence of strong stimulation. Figure 12(e), (f) show that for the chosen set of parameters 25 % of regulatory cells are enough to cancel out the high levels of IL2 as no activated cells (marked with gold) are found.
We observed that for a large secreting rate the surrounding cells are activated as well for the spatial range over which the IL2 signal can occur.
Result 16.
Secreting cells produce shortrange signals despite of fast diffusion. For the chosen set of parameters, the radius of activation covers 2 cell distances (20 μm). For the higher secretion rate in the Result 12 it covers 4 cell distances (50 μm). Only clusters of secreting cells will cause longrange signals [13].
Conclusion
In this work we presented the simulation results of new mathematical models for signaling pathways. Including the aspect of the diffusion of the molecules and the benefit of measured data, our models show a special structure described by a system of coupled PDE/ODE with Robin type boundary conditions. The three dimensional simulations in several realistic cell geometries with specific parameters show the dominance of the chemical reactions over diffusion in all models. The chemical reactions determine the gradient formation and their result to the cell fate and finally, the spatial range of the signal. Nevertheless, the diffusion must be considered in cells with irregular form or for slow diffusing molecules (D<1 μ m ^{2}/s).
Declarations
Acknowledgments
Part of this work was supported by the Helmholtz Alliance on Systems Biology (SB Cancer, Submodule V.7).
I acknowledge the financial support by Deutsche Forschungsgemeinschaft and RuprechtKarlsUniversität Heidelberg within the funding programme Open Access Publishing.
I would like to thank my advisor, Rolf Rannacher, who supervised the work and provided some of his students for this research topic, Rebecca Neumann, Juliane Claus and Daniel Gerecht for their numerical simulations. I thank Susanne Krömker, Lisa Kolb and Marcus Schaber for their contribution to the three dimensional visualization techniques. The results presented here would not have been possible without the contribution of my collaborators in experimental and theoretical Systems Biology: Kevin Thurley, Andrea C. Pfeifer, Thomas Höfer and Ursula Klingmüller. I would like to thank Rebecca Wade from the Heidelberg Institute of Theoretical Studies for the opportunity to communicate this subject at the Biological Diffusion and Brownian Dynamics Brainstorms. This paper was written within the framework of these meetings. Finally, I would like to thank Sara Lee for proof reading and her helpful comments.
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
 Kestler HA, Wawra C, Kracher B, Kühl M. Network modeling of signal transduction: establishing the global view. BioEssays. 2008; 30:1110–25.View ArticleGoogle Scholar
 Turing MA. The chemical basis of morphogenesis. Phil Trans R Soc B. 1952; 237:37–72.View ArticleADSGoogle Scholar
 Amonlirdviman K, Khare NA, Tree DRP, Chen WS, Axelrod JD, Tomlin CJ. Mathematical modeling of planar cell polarity to understand domineering nonautonomy. Science. 2005; 307:423–6.View ArticleADSGoogle Scholar
 Chaplain MAJ, Gerisch A. Mathematical modelling of cancer cell invasion of tissue: Local and nonlocal models and the effect of adhesion. J Theo Bio. 2008; 250(4):684–704.View ArticleMathSciNetGoogle Scholar
 Grimbs S, Arnold A, Koseska A, Kurths J, Selbig J, Nikoloski Z. Spatiotemporal dynamics of the calvin cycle: multistationarity and symmetry breaking instabilities. Biosystems. 2011; 103:212–23.View ArticleGoogle Scholar
 Brown GC, Kholodenko BN. Spatial gradients of cellular phosphoproteins. FEBS Lett. 1999; 457:452–4.View ArticleGoogle Scholar
 Caudron M, Bunt G, Bastiaens P, Karsenti E. Spatial coordination of spindle assembly by chromosomemediated signaling gradients. Science. 2005; 309(5739):1373–6.View ArticleADSGoogle Scholar
 Crank J. Mathematics of Diffusion. London: Oxford & Clarendon Press; 1975.Google Scholar
 Friedmann E, Pfeifer AC, Neumann R, Klingmüller U, Rannacher R. Interaction between experiment, modeling and simulation of spatial aspects in the jak2/stat5 signaling pathway In: Bock H, Carraro T, Jäger W, Körkel S, Rannacher R, Schlöder JP, editors. Model Based Parameter Estimation  Theorie and Application. Berlin Heidelberg: Springer: 2013. p. 125–43.Google Scholar
 Klingmüller U, Bauer A, Bohl S, Nickel PJ, Breitkopf K, Dooley S, et al. Primary mouse hepatocytes for systems biology approaches: a standardized in vitro system for modelling of signal transduction pathways. Syst Biol. 2006; 153:433–7.View ArticleGoogle Scholar
 Friedmann E, Neumann R, Rannacher R. Well posedness of a linear spatio temporal model of the jak2/stat5 signaling pathway. Comm Math Anal. 2013; 5(2):76–102.MathSciNetGoogle Scholar
 Claus J, Friedmann E, Klingmüller U, Szekeres T. Spatial aspects in the SMAD signaling pathway. J Math Biol. 2013; 67:1171–97.View ArticleMathSciNetMATHGoogle Scholar
 Thurley K, Gerecht D, Friedmann E, Höfer T. Threedimensional gradients of cytokine signaling between t cells. PLOS Comput Biol. 2015; doi:10.1371:1–22.Google Scholar
 Busse D, de la Rosa M, Hobiger K, Thurley K, Flossdorf M, Scheffold A, Höfer T. Competing feedback loops shape IL2 signaling between helper and regulatory T lymphocytes in cellular microenvironments. Proc Nat Acad Sci USA. 2010; 107:3058–63.View ArticleADSGoogle Scholar
 de la Rosa M, Rutz S, Dorninger H, Scheffold A. Interleukin2 is essential for CD4+CD25+regulatory T cell function. Eur J Immunol. 2004; 34:2480–8.View ArticleGoogle Scholar
 Scheffold A, Murphy KM, Höfer T. Competition for cytokines: T(reg) cells take all. Nat Immunol. 2007; 8:1285–7.View ArticleGoogle Scholar
 Thurley K. Numerische und analytische Lösungen eines räumzeitlichen Modells der Interleukin2 Expression in TLymphozyten. HumboldtUniversität zu Berlin: PhD thesis; 2007.Google Scholar
 GASCOIGNE. High Performance Adaptive Finite Element Toolkit. http://www.gascoigne.unihd.de. Accessed 2010.
 Bangerth W, Hartmann R, Kanschat G. Deal.II  a genera purpose object oriented finite element library. ACM Trans Math Softw. 2007; 33:24–12427.View ArticleMathSciNetGoogle Scholar
 Bangerth W, Rannacher R. Adaptive Finite Element Methods for Differential Equations. Basel: Birkhäuser; 2003.View ArticleMATHGoogle Scholar
 Becker R, Rannacher R. An optimal control approach to a posteriori error estimation in finite element methods. Acta Numerica. 2001; 10:1–102.View ArticleMathSciNetMATHGoogle Scholar
 Kolb L. Visualizing HighResolution Numerical Data with Isosurfaces using Topological Methods. University Heidelberg, IWR, Numerical Geometry: PhD thesis; 2013.Google Scholar
 Deheuninck J, Luo K. Ski and SnoN, potent negative regulators of TGFbeta signaling. Cell Res. 2009; 19(1):47–57.View ArticleGoogle Scholar
 Wilkinson DS, Ogden SK, Stratton SA, L PJ, Nguyen TT, Smulian GA, et al. A direct intersection between p53 and transforming growth factor beta pathways targets chromatin modification and transcription repression of the alphafetoprotein gene. Mol Cell Biol. 2005; 25(3):1200–12.View ArticleGoogle Scholar
 Huse M, Lillemeier BF, Kuhns MS, Chen DS, Davis MM. T cells use two directionally distinct pathways for cytokine secretion. Nat Immunol. 2006; 7:247–55.View ArticleGoogle Scholar
 PeronaWright G, Mohrs M. Sustained signaling by canonical helper T cell cytokines throughout the reactive lymph node. Nat Immunol. 2010; 11:520–6.View ArticleGoogle Scholar
 Feinerman O, Jentsch G, Tkach K, Coward J, Hathorn M, Sneddon MW, et al. Singlecell quantification of il2 response by effector and regulatory t cells reveals critical plasticity in immune response. Mol Syst Biol. 2010; 6:437.View ArticleGoogle Scholar
 Sabatos C, Doh J, Chakravarti S, Friedman R, Pandurangi P, Tooley AJ, et al. A synaptic basis for paracrine interleukin2 signaling during homotypic t cell interaction. Immunity. 2008; 29:6136–43.View ArticleGoogle Scholar