 Research Article
 Open Access
 Published:
PDE/ODE modeling and simulation to determine the role of diffusion in longterm and range cellular signaling
BMC Biophysics volume 8, Article number: 10 (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.
Background
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
Let x∈Ω and t∈ [0,T], 0≤T<∞ be two independent variables (spatial and temporal respectively) and u _{ i }(t,x) denote the spatialtemporal molecule concentrations (i=0,1,…,n) involved in the specific pathway at position x and time t, \(n\in \mathbb {N}\). Also let \(\,\Omega \subset {\mathbb {R}}^{3}\,\) be the space resolved part (e.g. cytoplasm or extracellular space) and thus bounded by a sufficiently smooth boundary ∂ Ω. We use the standard notation of spacetime function spaces. For a real function Banach space X with norm ∥·∥ on Ω, the space L ^{p}(0,T;X) consists of all measurable functions u:[0,T]→X with
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 general structure of our models looks like a coupled twocompartment model. In Ω we have reactiondiffusion equations (PDE) for each diffusing concentration u _{ i } with diffusion rate D _{ i }≥0, i=0,…,m :
The functional f describes the kinetical part of the law of mass action and the coupling with involved concentrations u _{ j } entering the domain Ω. We potentially have a nonlinear Robin boundary condition on a partial boundary Γ⊂∂ Ω :
For the mathematical correct formulation we impose an initial condition from measurement:
Under the assumption of the wellmixed molecule concentrations over a certain domain (for example in the nucleus in our intracellular signaling models or in the entire cell in our intercellular signaling models), we have ODE:
with the spaceindependent averaged value \(\bar {u}_{i}=\frac {1}{\partial \Omega }\int _{\partial \Omega } u_{i}(t,s)ds\) and initial conditions
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
We implement the coupled PDE/ODE system in the FiniteElementplatforms developed in our group (Gascoigne [18] and deal.II [19]). The principal components of our numerical methods are:

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
Our model of the JAK2/STAT5 signaling pathway and the parameters can be found in [9]. The STAT5 molecules form dimers in both activated and nonactivated states. If our concentration variables describe the dimers instead of the monomers, the resulting model equations are linear with a unique solution [11]. We consider the model for different cell types and different cell geometries:

a mathematical model of fibroblast (NIH3T3, Figs. 1, 3 and 6)

a reconstructed fibroblast via microscopy data (Figs. 1 and 4)

an artificial cell where we added two cytoplasmic extensions (filopodia) to the reconstructed geometry (Fig. 5).
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.
To examine the spatial distribution of the diffusive molecules, activated and nonactivated STAT5 in the cytoplasm (u _{0},u _{1}), we evaluate the maximum relative variation of a concentration at a given time point t:
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.
The results of the simulation of the JAK2/STAT5 model are presented in the Figs. 2, 3, 4, 5 and 6.
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.
A cell has a complex structure thus the diffusion of the molecules therein is not isotropic like in homogeneous materials, i.e. the same diffusion in every direction. In some parts of the cytoplasm there is a retardation in the motion of molecules possibly caused by binding processes and a speedup by active transport or aided diffusion. Especially in signaling pathways, the molecules are recruited towards the nucleus. Often they are transported by motorproteins along the microtubules, so that the diffusion can be faster in direction to the nucleus. To model this biological behavior we consider the convectiondiffusion equation with a time and space dependent diffusion tensor \(\tilde {D}(t,x)\). For i=0,…,m and j=m+1,…,n, we have
v is the velocity of the motorproteins. To determine v it requires additional experiments. As a first approximation, we use the anisotropic diffusion equations with an inhomogeneous but constant diffusion tensor for each diffusing species:
with
For the isotropic diffusion in Eq. (2) we used the same diffusion coefficient for the activated and nonactivated STAT5, \(D_{i}=\bar {D}\), i=0,1.
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.
As a second intracellular signaling pathways we analyze the SMAD signaling pathway in hepatocytes. Details of the model and the biological function of the pathway can be found in [12]. The SMADmolecules undergo a more complex oligomerization. They form dimers and trimers and react with other molecules building multicomplexes that are described by nonlinear and more complex equations. For the cell shape we also consider various geometries:

a mathematical model of hepatocyte (Fig. 7 left)

a reconstructed geometry from microscopy data (Fig. 7 central)

an artificial geometry assembled with two cytoplasmic extensions (Fig. 7 right).
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].
To observe the influence of the cytosolic gradients on any further processes concerning gene regulation we analyze the difference in the nuclear trimer concentration which binds to the DNA. Figure 8 shows the result as follows:
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.
For the further investigation on this particular case where the greater gradients in the cytoplasm have only small effect on the amount of concentration in the nucleus we compare the total amount of SMADtrimers in the two compartments (Fig. 9):
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.
A slower diffusion of the cytosolic trimer with D=1 μ m ^{2}/s instead of D=15 μ m ^{2}/s implied a 70 % cytosolic SMADtrimer concentration variation in the steady state (Fig. 10) compared to 18 % as seen in the previous result.
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].
For numerical calculations we choose a part of the three dimensional lymph node with 125 or 218 cells (Fig. 11). 25 % of these cells are randomly chosen as secreting IL2 (secretory T cells) and the rest of the cells are absorbing IL2 (T helper and regulatory cells). The cells compete for IL2 and those who absorb more will upregulate their receptors and consume even more IL2. So called regulatory T cells have a higher absorbing rate and may downregulate the signal, for example to avoid autoimmune reactions. The interaction of the cells generates a space and time dependent dynamics which describes the competition of the cells for IL2. This dynamics depends on the position of the different cell types and stabilizes after 30 h in an inhomogeneous steady state. The winner cells are activated as a large number of receptors (more than 4000 receptors after 30 h of simulation time) are formed on the surface.
Result 10.
In Fig. 12 the particular isosurfaces of the IL2 distribution are presented in the extracellular space at specific time points (after 1, 9 and 30 h of activation) obtained by topological methods [22]. We confirmed the fastdiffusivity of IL2; after only an hour IL2 diffused everywhere in a high concentration. Then the local chemical reactions on each cell start to determine the behavior of the dynamics. If a cell has a greater chance to absorb more IL2 it express more receptors. Due to the higher amount of expressed receptors the cell absorbs even more IL2 so that in its vicinity there is not enough IL2 left for the other cells to be activated. Consequently, an inhomogeneous IL2 pattern and high concentration gradients are formed. In our model we have 31 secretory cells out of 125. After 9 h of activation we observed 20 activated cells and 33 cells after 30 h in the steady state. The same simulations with 31 regulatory T cells showed an IL2 concentration of one order of magnitude (ten times) lower after one hour (comparison Fig. 12 a) with d)). After 9 h we observed a very low IL2 concentration and no activation. This behavior remains until the steady state. We can conclude that in our model 31 regulatory T cells can completely downregulate the signal.
Result 11.
The Fig. 13 left shows the IL2 distribution on the surface of 218 cells. The concentration varies over three orders of magnitude. The Fig. 13 right shows the expression level of the cells: 60 out of 218 cells have more than 4000 IL2 receptors and are thus considered as activated. By introducing 54 regulatory cells into the lymph node this activation would be suppressed completely.
To analyze the range of action of the cytokine IL2 we constructed a three dimensional in silico experiment containing 3375 cells (Fig. 14). The cell in the middle is the single secreting cell in the system and all others are naive T helper cells which can be activated when enough IL2 is available.
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.
Figure 10 demonstrates that different sizes of gradients in the cytoplasm (70 %, 18 % respectively) caused by two different diffusion coefficients (D=1 μ m ^{2}/s, D=15 μ m ^{2}/s) give rise to the same dynamics in the nuclear trimer concentration. For the considered set of parameters our numerical simulations show that the diffusion of molecules in the cytoplasm always plays a role in cells with an irregular shape like dendritic cells or cells with cytoplasmic extensions. In case the diffusion coefficient D is of a few orders of magnitude lower than in the considered pathways (D<<15 μ m ^{2}/s), visible effects may appear in the dynamics also for cells with regular shape. The deceleration of diffusion of a factor 10 does not have an effect (Figs. 10 and 15). This gives an idea that for our applications the measurement of diffusion coefficients allows more flexibility.
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.
Figure 15 shows that the cell with extensions exhibits a greater relative concentration variation in the cytoplasmic concentration (45 %) and also a visible deviation in the nuclear trimer concentration. The molecules are activated on the outer membrane and need a longer time to travel through the cytoplasmic dendrite to the nucleus. This implies that less molecules are available for the processes in the nucleus, and diffusion can change gene expression in a cell with irregular shape. This is explained by the inhibition of pSMAD2 and trimer nuclear accumulation in [12]. SMAD2 molecules persist longer in the cytoplasmic domain. Therefore, the chance for other proteins to compete for SMAD2 binding increases. In the end less molecules will be available in the nucleus to bind to the DNA. Nevertheless, SMAD molecules are latent transcription factors that interact with a plethora of different coactivators, corepressors and polymerase. One of the wellstudied cofactors that binds to nuclear SMAD trimers is SnoN. Recent data suggests that SnoN modulates the effect of TGF beta [23]. In Fig. 16 we compare the SnoN and the multiprotein complex SnoNSMADtetramer concentration from the ODE with the PDE model in the cell geometry with extensions. Due to the cytosolic trimer diffusion, there are less trimers available in the nucleus than for the ODE model. So, in the nucleus more SnoN molecules and less complexes than in the ODE model can develop in time. The deviation in the complex concentration is present only from 15 to 50 min after activation which implies that due to nuclear accumulation there are enough trimer molecules available in the nucleus for complex formation when enough time is provided. The amount of the concentration deviation seems small: with the values of 0.1 nM (2 % of the total concentration) for SnoN and 0.05 nM (6 % of the total concentration) for the complex SnoNSMADtetramer. Nevertheless, due to the feedback loops even small changes in nuclear SMAD trimers can have a profound effect on the multifactor complexes. SnoN interacts also with the transcription factor p53 [24]. In conclusion the temporal availability of SMADtrimers due to the modulatory effect of coactivators and corepressors can result in altered cell fate in cells with irregular shape.
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).
References
 1
Kestler HA, Wawra C, Kracher B, Kühl M. Network modeling of signal transduction: establishing the global view. BioEssays. 2008; 30:1110–25.
 2
Turing MA. The chemical basis of morphogenesis. Phil Trans R Soc B. 1952; 237:37–72.
 3
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.
 4
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.
 5
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.
 6
Brown GC, Kholodenko BN. Spatial gradients of cellular phosphoproteins. FEBS Lett. 1999; 457:452–4.
 7
Caudron M, Bunt G, Bastiaens P, Karsenti E. Spatial coordination of spindle assembly by chromosomemediated signaling gradients. Science. 2005; 309(5739):1373–6.
 8
Crank J. Mathematics of Diffusion. London: Oxford & Clarendon Press; 1975.
 9
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.
 10
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.
 11
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.
 12
Claus J, Friedmann E, Klingmüller U, Szekeres T. Spatial aspects in the SMAD signaling pathway. J Math Biol. 2013; 67:1171–97.
 13
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.
 14
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.
 15
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.
 16
Scheffold A, Murphy KM, Höfer T. Competition for cytokines: T(reg) cells take all. Nat Immunol. 2007; 8:1285–7.
 17
Thurley K. Numerische und analytische Lösungen eines räumzeitlichen Modells der Interleukin2 Expression in TLymphozyten. HumboldtUniversität zu Berlin: PhD thesis; 2007.
 18
GASCOIGNE. High Performance Adaptive Finite Element Toolkit. http://www.gascoigne.unihd.de. Accessed 2010.
 19
Bangerth W, Hartmann R, Kanschat G. Deal.II  a genera purpose object oriented finite element library. ACM Trans Math Softw. 2007; 33:24–12427.
 20
Bangerth W, Rannacher R. Adaptive Finite Element Methods for Differential Equations. Basel: Birkhäuser; 2003.
 21
Becker R, Rannacher R. An optimal control approach to a posteriori error estimation in finite element methods. Acta Numerica. 2001; 10:1–102.
 22
Kolb L. Visualizing HighResolution Numerical Data with Isosurfaces using Topological Methods. University Heidelberg, IWR, Numerical Geometry: PhD thesis; 2013.
 23
Deheuninck J, Luo K. Ski and SnoN, potent negative regulators of TGFbeta signaling. Cell Res. 2009; 19(1):47–57.
 24
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.
 25
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.
 26
PeronaWright G, Mohrs M. Sustained signaling by canonical helper T cell cytokines throughout the reactive lymph node. Nat Immunol. 2010; 11:520–6.
 27
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.
 28
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.
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.
Author information
Additional information
Competing interests
The author declares that he has no competing interests.
Authors’ contributions
The author declares to be the principle investigator of the work presented here. The idea to analyze the role of diffusion in cellular signaling was worked out during my postdoc time in the Numerical Simulation group of Rolf Rannacher. The concept of the work was developed due to my interdisciplinary education in applied mathematics at the Interdiciplinary Center of Scientific Computing in Heidelberg. Some of the numerical computations were done under my supervision in the frame of diploma theses and one dissertation which should be finished soon.
Rights and permissions
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.
About this article
Cite this article
Friedmann, E. PDE/ODE modeling and simulation to determine the role of diffusion in longterm and range cellular signaling. BMC Biophys 8, 10 (2015) doi:10.1186/s1362801500248
Received
Accepted
Published
DOI
Keywords
 Systems of coupled differential equations
 Reactiondiffusion systems
 Numerical simulation
 Cellular signaling