Spatial modeling of the membrane-cytosolic interface in protein kinase signal transduction

The spatial architecture of signaling pathways and the interaction with cell size and morphology are complex, but little understood. With the advances of single cell imaging and single cell biology, it becomes crucial to understand intracellular processes in time and space. Activation of cell surface receptors often triggers a signaling cascade including the activation of membrane-attached and cytosolic signaling components, which eventually transmit the signal to the cell nucleus. Signaling proteins can form steep gradients in the cytosol, which cause strong cell size dependence. We show that the kinetics at the membrane-cytosolic interface and the ratio of cell membrane area to the enclosed cytosolic volume change the behavior of signaling cascades significantly. We suggest an estimate of average concentration for arbitrary cell shapes depending on the cell volume and cell surface area. The normalized variance, known from image analysis, is suggested as an alternative measure to quantify the deviation from the average concentration. A mathematical analysis of signal transduction in time and space is presented, providing analytical solutions for different spatial arrangements of linear signaling cascades. Quantification of signaling time scales reveals that signal propagation is faster at the membrane than at the nucleus, while this time difference decreases with the number of signaling components in the cytosol. Our investigations are complemented by numerical simulations of non-linear cascades with feedback and asymmetric cell shapes. We conclude that intracellular signal propagation is highly dependent on cell geometry and, thereby, conveys information on cell size and shape to the nucleus.

Introduction Cells need to respond to a large variety of external stimuli such as environmental changes or extracellular communication signals. Signals transmitted from cell surface receptors to target genes in the nucleus are frequently transduced by cascades of covalent protein modifications. These modifications consist of inter-convertible protein forms, for instance, a phosphorylated and an unphosphorylated protein. Signaling cascades occur in many different variations including mitogen-activated protein-kinase (MAPK) cascades and small GTPase cascades.
Signal transduction mechanisms carried out by networks of protein-protein interactions are highly modular and regulatory behavior arises from relatively simple modifications [1]. The spatial arrangement of signaling cascades varies in different biological systems. We focus on the localization of signaling components, which can be tethered to the cell-membrane or freely diffuse in the cytosol. Tethering of signaling molecules to the cell-membrane can be mediated by lipidation modifications [2][3][4][5][6], co-localization by membrane-bound scaffolds [7] or membrane anchoring proteins [8]. Frequently, the first steps of signal transduction occur at the membrane and are then continued into the cytosol. We investigate linear signaling cascades with different realizations of spatial arrangements of signaling components as shown in Fig 1. Here, we focus on the membrane-cytosolic interface, which is included in the signaling motif shown in Fig 1(B) and 1(C).
In many experimental and theoretical studies on signaling cascades, the cell is regarded as a number of well-mixed compartments with no variation in size, shape or organelle location. Attempts of a quantitative description of signaling cascades with a focus on temporal aspects have been made in [9][10][11][12]. However, the spatial description of signaling processes has received less attention despite its relevance in understanding cell morphology and growth regulation in time and space [13]. Examples of spatial effects on the length scale of single cells range from the yeast mating process [14,15] to the propagation of spatial information in hippocampal neurons which is controlled by cell shape and vice versa [16,17].
Since the cytosol scales with cell volume and the cell membrane with the cell surface area, reactions on the membrane and in the cytosol scale with the cell-surface to cell-volume ratio. For instance, we obtain an area/volume ratio of / 3/R cell for a spherical cell geometry, where R cell is the cell radius. We will show that this scaling affects the global phosphorylation rate of signaling proteins that diffuse in the cytoplasmic volume, which depends on cell size. While cytosolic gradients naturally occur from the membrane to the nucleus, membrane-bound components can only form gradients along the membrane, which changes the response to heterogeneous signals. Furthermore, the diffusion on the membrane is much slower for membrane-bound components than for cytosolic components [18]. Both of these factors are expected to largely change signal transduction properties of the pathway.
An analysis and comparison of spatial signal transduction motifs in response to spatially homogeneous and heterogeneous signals is presented in this study. The natural extension of widespread used ordinary differential equations are bulk-surface partial differential equations [19,20]. Here, bulk refers to the cellular compartments that are represented as a volume such as the cytoplasm or the nucleus, while surface refers to all cellular structures that are represented as an area such as the cellular or nuclear membrane. Since their introduction to cell signaling systems [21], bulk-surface partial differential equations have been successfully employed in several models for cell polarization [18,[22][23][24]. However, membrane-cytosolic interfaces at different stages of a signaling cascade have not yet been investigated.
We start with an analysis of two different motifs with simplified linear kinetics, which allows to develop exact analytical solutions of the steady state. Both motifs differ in their cell size dependence and we show further that their behavior can be drastically different from the assumption of well-mixed compartments. The time-scaling of signal transduction is investigated using the method of local accumulation times [25]. We continue by investigating the response and sensitivity to spatially heterogeneous signals such as signaling gradients for symmetrical and asymmetrical cell shapes. In the last section, we proceed with numerical investigations of systems with negative feedbacks which lead to cell-size dependent oscillations. A Fourier analysis in time is used to provide insight into the dependency of oscillation frequency and amplitude on cell size. Depending on the spatial motif, cell size limits for the extinction of oscillatory behavior are obtained.
We start with a linear signaling cascade with different localizations of the membranecytosolic interface as shown in Fig 1. We employ a simple cascade model from [9], in which stimulation of a receptor leads to the consecutive activation of several down-stream protein kinases. This model is extended into space in the following. We assume a linear cascade with N components, where the first M < N components are localized at the membrane while the remaining N − M components are assumed to freely diffuse in the cytosol. The equations for the membrane-bound components read @P n @t ¼ D mem D G P n þ v a n À v d n on the cell membrane; for n ¼ 1; . . . ; M: Here, P 1 ðx; tÞ; . . . ; P M ðx; tÞ are the local concentrations of signaling molecules on the cell membrane. The activation rate of the first signaling component v a 1 is assumed to be dependent on the input signal, which is denoted by P 0 ðx; tÞ. The input signal on the cell surface can be a trigger on the cell membrane or arise from an extracellular signal. All of these species are functions of space and time, wherex is a point on the membrane and t is the time. Diffusion along the cell membrane, which is assumed to be a two-dimensional curved surface in three-dimensional space, is described by the Laplace-Beltrami operator Δ Γ and the diffusion coefficient D mem . Since the membrane is a surface in three-dimensional space with negligible thickness, the natural unit for concentrations of the cell membrane-bound species P n (n = 1, . . ., M), is molecules per area. Molecular concentrations of signaling molecules are frequently provided in nanomolar or micromolar (nM or μM). For convenience, we therefore use the units nanomolar or micromolar times micrometer (nMμm or μMμm) for the membrane-bound signaling molecules. Note that 1 μMμm % 602 molec/μm 2 .
The phosphorylation rates v a n as well as the dephosphorylation rates v d n have units molecules per area and time. If the input signal is homogeneous in space, meaning P 0 ðx; tÞ ¼ P 0 ðtÞ, all spatial fluxes D mem r Γ P n are zero and the equation system for the membrane-bound species can be described by an equivalent system of ordinary differential equations (S1 Appendix). In contrast to the membrane-bound signaling components P 1 , . . ., P M , the signaling component P M+1 can freely diffuse in the cytosol. For the modeling of the membrane-cytosolic interface, we need to include diffusion in the cytosol and reactions on its boundaries, which are the membranes. These processes are modeled by a reaction-diffusion equation with the boundary condition Since P M+1 is activated by the upstream component P M , which is tethered to the membrane, there is a phosphorylation reaction only at the cell membrane but not in the interior of the cytosol. This reaction is, therefore, modeled as a boundary condition. The reactions at the membrane-cytosolic interface are described by the phosphorylation rate v a Mþ1 and the inactivation rate v i , both with units molecules per area and time. The species P M+1 diffuses freely in the cytosolic volume with the diffusion rate D cyt and therefore its local concentration is described in units molecules per volume. The dephosphorylation rate v d Mþ1 in the cytosol is given in molecules per volume and time. Note that the inactivation rate and v i can be invoked by membrane-bound phosphatases or saturation of phosphorylation at the membrane. Both, v a and v i , comprise the kinetics at the membrane-cytosolic interface. For the flux on all other membrane enclosed organelles we assume a zero-flux condition The equations for the components of the downstream cytosolic cascade read @P n @t ¼ D cyt DP n þ v a n À v d n ; in the cytosol; for n ¼ M þ 2; . . . ; N: The concentrations of the cytosolic components at positionx in the cytosolic volume at time t are described by functions P n ðx; tÞ with n = M + 1, . . ., N. For the cytosolic components we assume zero-flux conditions: À D cyt rP n Áñ ¼ 0; on the cell membrane; ð6Þ À D cyt rP n Áñ ¼ 0; on the nuclear membrane; In classical MAPK cascades the last component of the cascade, which is the phosphorylated MAPK, is imported into the nucleus. Examples range from Hog1 nuclear import in yeast [26,27] to the import of ERK in mammals [28]. In this case, the boundary condition Eq (7) on the nucleus for the last cytosolic component P N needs to be modified to where represents a nuclear-import reaction rate on the nuclear membrane. Unless otherwise stated, a zero-flux boundary condition is assumed on the nucleus throughout this paper. We will test and compare systems with three components N = 3 as shown in Fig 1, where the spatial arrangement of the components is varied. Here, M = 2 describes the case of two membrane-bound and one cytosolic element (motif Fig 1B) and M = 0 the case of only cytosolic components (motif Fig 1C). In the following the case M = 2 is referred to as mixed membrane-cytosolic (MMC) and M = 0 as pure cytosolic (PC) cascade.

The mixed membrane-cytosolic cascade is strongly size dependent
We start this section with an analysis of a spherical cell and then generalize the analysis to arbitrary cell shapes. A spherical cell of radius R cell with a spherical nucleus of radius R nuc placed in the center of the cell is assumed in the following.
The input signal is denoted by P 0 (t) and is assumed to be homogeneous on the cell surface. The concentrations of protein kinases are described by functions P i (r, t) depending on space and time. Note, since the cellular geometry is radially symmetric and the input signal P 0 acts homogeneously on the cell membrane, these functions depend only on the radial distance from the cell center, denoted by r, and time t. In the following analysis, the kinetic rates are linearized, meaning that we assume v a n ¼ a n P nÀ 1 and v d n ¼ b n P n for the phosphorylation and dephosphorylation, respectively. The inactivation rate v i at the membrane-cytosolic interface is as well linearized by v i = γP M+1 . Note that the interface kinetics can be reformulated as , from which it can easily be seen that the activation at the membrane saturates at P Mþ1 ¼ a Mþ1 g P M . The model equations for the mixed membrane-cytosolic cascade (MMC) with linearized kinetics read and boundary conditions for the cytosolic species P 3 : À D cyt rP 3 Áñ ¼ a 3 P 2 À gP 3 on the membrane; ð12Þ À D cyt rP 3 Áñ ¼ À P 3 at the nucleus: ð13Þ There are several estimates of phosphatase activity and diffusion coefficients for MAPK signaling components. The diffusion coefficient, D cyt , of globular cytosolic proteins has been shown to be in the range 1 -10 μm 2 s −1 , while the diffusion coefficient of membrane-bound components, D mem , is much lower with a value in the range of 10 −3 -0.1 μm 2 s −1 [29][30][31]. The phosphatase rates β n range over three orders of magnitude 0.1 -100 s −1 [32,33]. In the case of Fus3, which is the MAPK in the mating pathway of the yeast S. cerevisiae, the diffusion coefficient and cytosolic dephosphorylation rate were estimated to be 4.2 μm 2 s −1 and 1 s −1 , respectively [14]. See Table 1 for an overview on parameter values and units.
We begin with a steady state analysis of this system in the parameter regimes of interest and assume that the signal P 0 is constant over time. Here and in the following we indicate the steady distribution with a bar, meaning that " P n denotes the steady state of P n . The steady state of the first two elements is given by " For the steady state of P 3 , the solution is given by where i 0 and k 0 are modified spherical Bessel functions of the first and second kind, respectively [37]. Note that i 0 is increasing with r (distance from the cell center), while k 0 is a decreasing function of r. The coefficients A and B are derived in S1 Appendix. If we neglect the nucleus or there is no nucleus in the cytosol, meaning R nuc = 0, the coefficient B becomes zero. The steady state solution for different cell sizes is shown in Fig 2. The concentration is maximal at the cell membrane and decays towards the nucleus. An estimate of the decay length L gradient of the intracellular gradient (with highest concentration at the membrane) is given by . This decay length can be compared with the actual cell size. Their ratio is called the Thiele modulus, a dimensionless measure defined [33]. For F ) 1 strong intracellular gradients and concentration heterogeneities of signaling molecules are to be expected, while for F ( 1 the concentration is almost homogeneous. Since the Thiele modulus relates the diffusion coefficient and degradation rate to cell size, it is an important parameter to investigate gradient formation [33] and signal propagation for several cascade levels [38]. However, in a three-dimensional space L gradient can not be interpreted as the actual gradient anymore, since its derivation is Spatial modeling of the membrane-cytosolic interface in protein kinase signal transduction based on the assumption of a one-dimensional geometry. In addition, if an excluding volume such as the nucleus is assumed, the one-dimensional L gradient overestimates the concentration gradient. For example, if we assume D cyt = 4.0 μm 2 /s and β = 1.0 s −1 , we obtain L gradient = 2.0 μm. In the case of the classical one-dimensional simplification a decay proportional to / exp(−x/L gradient ) is assumed, which suggests a concentration decrease in a distance of x = 2 μm by a factor of exp(−x/L gradient ) % 0.37. However, in a spherical cell with radius R cell = 3 μm with excluding volume R nuc = 1 μm, the concentration decreases only by a factor of 0.77 in a distance of 2 μm from the cell membrane. The effect of cell size on intracellular concentration gradients is shown in Fig 2. The cell size dependence in cell signaling systems does not only arise by the characteristic length scale for intracellular gradient formation, but by the change of average intracellular concentration levels with cell size. We start with the simplifying assumption that there is no nucleus or excluding volume in the cytosol, meaning R nuc = 0. In this case the steady state solution reads The modified spherical Bessel functions i 0 and i 1 are monotonically increasing functions with lim for cells with small F and, therefore, the phosphorylation reaction at the membrane is at saturation in this case. For large F, we obtain from lim These estimates also hold in the case of an excluding volume which is the nucleus and we obtain the estimate for the concentration " P 3 ðR cell Þ at the cell membrane: The dependence of absolute concentration levels on the membrane-cytosolic interface is shown in Figs 2 and 3 for a set of different parameters. For a large inactivation rate at the mem- , the cell size dependence decreases. Therefore, cell size dependence is mainly determined by γ and ffiffiffiffiffiffiffiffiffiffiffi ffi D cyt b 3 q but is independent of the phosphorylation rate α.
We can further investigate the evolution of the average concentration levels, which depends on the concentration at the cell membrane and the strength of the intracellular gradient. In case of arbitrary cell shapes with cell volume V cell and cell membrane area M cell , the average concentration is obtained from In [33], analytical solutions for the average concentration in a spherical cell and a slab have been derived as functions of the Thiele modulus. However, since the derivations are restricted to cell geometries, where an explicit analytical solution of the reaction diffusion equation is  available, we introduce an alternative approach to estimate average concentration levels depending on the cytosolic volume and cell membrane area. Since the signal propagates from the cell membrane to the cytosol, the cell membrane can be regarded as a source, while the cytosolic volume, where phosphorylated signaling molecules are dephosphorylated, can be regarded as a sink. This idea can be derived mathematically by integration of Eq (2) and application of Green's theorem, which results in Here, the production on the left hand side of the equation depends on the cell membrane area, which is balanced by the degradation in the cytosol on the right hand side of the equation. On the basis of the equation of mass conservation (see S1 Appendix) in reaction diffusion systems, we introduce the following measure: This measure has the property L Mþ1 ¼ " P avg Mþ1 for γ = 0 or β = 0, which holds for arbitrary cell shapes. Furthermore, for a spherical cell the estimate " P avg

Mþ1
L Mþ1 " P max Mþ1 holds for β > 0 and γ > 0 (see S1 Appendix and Fig 3). Therefore, we use Λ M+1 as a proxy for the average concentration for arbitrary cell shapes, which can be easily calculated.
A comparison of the estimate Λ 3 to the average concentration is shown in Fig 3 for different cell shapes, which are a spherical cell, a rod shaped cell, a cell with one protrusion and a cell with two protrusions. These cell shapes occur for example in S. cerevisiae, S. pombe and haploid S. cerevisiae stimulated with mating pheromone [39,40]. The MMC cascade was simulated for these shapes with varying cell size. The measure Λ 3 is an exact predictor for the average concentration in the case γ = 0 for all cell shapes and slightly overestimates the average concentration for γ = 1 μm s −1 .
Furthermore, we investigated the concentration differences of " P 3 between cell membrane and nucleus and compared them to the average concentration in the cytosol. For the spherical cell the average concentration levels of " P 3 in the cytosol as well as on the membrane were the lowest, which is expected since the surface to volume ratio is the lowest among all shapes. The concentration differences between membrane and nucleus were the highest for the spherical cell and the cell with two protrusions and the average concentration at the nucleus decreased almost to zero for large cells. For the rod shape cell the concentration differences were the smallest, since the distance along the short axis is small and the concentration does not drop as sharply as for the other cell shapes.
We furthermore established a correspondence to the evolution of the average concentration levels in time. In the case γ = 0 and for arbitrary cell shapes, the average concentration levels Spatial modeling of the membrane-cytosolic interface in protein kinase signal transduction follow the system of ordinary differential equations dP avg 1 ðtÞ dt ¼ a 1 P avg 0 ðtÞ À b 1 P avg 1 ðtÞ; ð20Þ where P avg 1 and P avg 2 are the average concentration levels in molecules per cell membrane area. This system of ordinary equations can be obtained by integrating Eqs (9)-(13) over their respective spatial domains. See S1 Appendix for details of the derivation. The steady state for the average concentration of " P 3 is given by Therefore, the average concentration level scales with the ratio of membrane area to cytosolic volume which is given by jM cell j jV cell j . The effective global phosphorylation rate for the average concentration of active signaling molecules in the cytosol is therefore determined byã 3 ¼ jM cell j jV cell j a 3 . These relations give us a correspondence between widespread used ordinary differential equations and the bulk-surface partial differential equations employed in this paper. In summary, we have strong cell size dependence, with decreasing concentrations for larger cells.

Efficient cytosolic transport via cytosolic cascades
In the following we consider a pure cytosolic (PC) cascade with three elements, in which all elements diffuse freely through the cytosol. The reaction-diffusion system is given by with boundary conditions on the membrane and at the nucleus Note that the membrane-cytosolic interface occurs at the first cascade level, meaning that only P 1 is activated at the membrane with rate α 1 P 0 − γP 1 . In the special case of β 1 = β 2 = β 3 = β analytical approximations to cytosolic cascades in a one-dimensional system have been derived in [41,42]. While a one-dimensional cellular geometry can be used to study gradient formation qualitatively, spatial effects such as the cell surface to volume ratio are neglected. Therefore, we derived exact analytical solutions to the linear system in three dimensions. The steady state solutions for " P n ðrÞ are expanded as follows The algebraic expressions of the coefficients A n, k and B n, k and their derivation are given in the S1 Appendix. In comparison to the MMC cascade, which was discussed in the previous section, the third cascade element P 3 is more evenly distributed in the cell and concentration gradients are much more shallow (see Fig 2). In order to quantify the concentration differences in a cell of arbitrary shape, we measure the concentration variance in the cell. Therefore, we introduce the variance as This measure has a close correspondence to the variance in image analysis [43,44]. In contrast to image analysis, where the square of the deviation of the fluorescence intensity from the average fluorescence intensity of a marker is integrated pixel-wise, the integration here is continuous. As for the analog in image analysis, the normalized variance is calculated as a measure for the deviation from the average [44,45]. While this measure is frequently used in auto-focus algorithms in image analysis, we suggest the normalized variance as a measure for the degree of localization of signaling molecules within a cell. An estimate the propagation of the normalized variance in the cytosolic cascade is given by Here, C s is a constant depending on cell shape and d is the cell diameter. Note that C = C s d is the Poincaré constant from the well known Poincaré inequality [46]. For convex cell shapes the estimate holds for C s ¼ 1 p . In the case of a convex cell shape of a small cell as yeast (without protrusion), we therefore have the estimate C n % 0.3 1 for D cyt = 3 μm 2 /s, α n = β n = 1 s −1 and a cell diameter of d = 6 μm. For this parameter set, the normalized variance decreases at least by 70% at the second cytosolic cascade level and by 90% at the third cytosolic cascade level (compared to the first cascade element). In general, for C n < 1, the normalized variance of the intracellular concentrations decreases with increasing cascade level and concentration differences in the cell are balanced out (see Fig 3).
Similar to the previous section, we derived an estimate for the average concentration. In this case, we employed the estimate Λ 1 to " P avg 1 , since P 1 is the cascade element that is activated at the membrane-cytosolic interface for the PC cascade. The average concentration " P avg 3 is related to " P avg 1 by " P avg 1 and, therefore, we used the approximation a 2 a 3 b 2 b 3 L 1 for " P avg 3 . As in the case of the MMC, this approximation is exact for γ = 0 and overestimates the average concentration sligthly for γ = 1 μm s −1 (see Fig 3).
Exact expressions for the steady states of the average concentration of signaling components in the case γ = 0 and for arbitrary cell shapes are given by Therefore, the average concentration of the third cascade element " P avg 3 takes the same values in the MMC and PC cascade. The major distinction of both spatial motifs is given by the fact that the concentration differences obtained at the cell membrane and nucleus are larger in the MMC cascade than in the PC cascade. Similarly as in the previous section, we can formulate a system of ordinary differential equations for the time evolution of average concentrations dP avg 1 ðtÞ dt ¼ jM cell j jV cell j a 1 P avg 0 ðtÞ À b 1 P avg 1 ðtÞ; ð35Þ dP avg 2 ðtÞ dt ¼ a 2 P avg 1 ðtÞ À b 2 P avg 2 ðtÞ; ð36Þ The dependence of absolute concentration levels on the membrane-cytosolic interface is shown in Fig 2 and figure in S1 Appendix for a set of different parameters.

The timing of spatial signaling
Time-resolved image-based analysis has shown that MAPK signaling pathways respond with a measurable signal in the nucleus in time scales of seconds to a few minutes. The Hog1 pathway response (phosphorylated MAPK) in budding yeast is at about 80% of its maximal activity within a minute [26]. Another example is the Src activation/deactivation cycle, where oscillations and pulses take place in the regime of seconds [47]. The timing of signal transduction in linear signaling cascades for well-stirred homogeneous systems has been analyzed in [9]. They concluded for weakly activated signaling cascades that phosphatases have a more pronounced effect than kinases on the rate and duration of signaling, whereas signal amplitude is controlled primarily by kinases. A thorough analysis of linear models assuming a homogeneous distribution of signaling molecules for different kinds of external stimuli has been recently worked out in [12]. We extended and compared these findings to spatial signal transduction omitting the simplification of homogeneous concentrations. How long does it take to establish an intracellular concentration gradient? How does diffusion change the timing of signal transmission from the cell membrane to nucleus? Which concentration differences are expected until a steady state is established? The time-scale analysis for spatial models is more difficult than for models based on a will-mixed assumptions due to high computational costs. Therefore, we used the recently introduced measure of accumulation times [25,48]. The approach of P n (r, t) to its steady state " P n ðrÞ at radial distance r from the cell center and time t can be characterized using the local relaxation function r n ðr; tÞ ¼ ð " P n ðrÞ À P n ðr; tÞÞ= " P n ðrÞ: The difference ρ n (r, t 1 ) − ρ n (r, t 2 ) can be interpreted as the fraction of the steady state level " P n ðrÞ that accumulated in the time interval [t 1 , t 2 ]. In an infinitesimal time interval [t, t + dt] the fraction of accumulated activated signaling molecules at steady state is given by À @r n ðr;tÞ @t dt.
The local accumulation time is defined as [25] t n ðrÞ ¼ À Z 1 0 t @r n ðr; tÞ @t dt: The accumulation time can be derived from the steady state solution even if no closed form of the time-dependent solution is known [25]. The timing of the average concentrations given in the system of ordinary differential equations for the MMC cascade (20)- (22) and the PC cascade (35)-(37) are the same and can be analytically expressed as This expression also coincides with signaling times calculated by Heinrich et al. [9]. However, for the spatial model the local accumulation times at the membrane and nucleus differ. The accumulation is generally faster at the membrane and slower at the nucleus, where the degree of the difference increases with cell size (see Fig 4). Furthermore, the two spatial motifs show significant differences. For the MMC cascade the accumulation time for the second element P 2 is exactly 1 b 1 þ 1 b 2 on the membrane, while it is shorter for the cytosolic species (compare Fig 4). The accumulation time of P 3 at the nucleus is, as expected, much longer. For small cells the intracellular concentration is spatially homogeneous and the approximation 1 b 1 þ 1 b 2 þ 1 b 3 holds, while the time for signal propagation to the nucleus increases with cell size. An analytical solution of the accumulation times for P 3 for the MMC cascade and the special case of R nuc = 0 can be derived [49], which is given in the S1 Appendix. However, for larger cells, the time for signal propagation to the nucleus increases with cell size. For the PC cascade, the increase in accumulation time at the nucleus with cell size is less pronounced than for the mixed-membrane cytosolic cascade.
While a constant stimulus was applied to calculate the accumulation times, we also tested a decaying signal P 0 ðtÞ ¼ P max 0 exp ðÀ ltÞ, with P max 0 ¼ 100 nMmm and solved the time-dependent system numerically. A comparison of the MMC and PC is shown in Fig 4. Interestingly, the concentration level at the membrane for the PC cascade decreases from the first cascade species P 1 to the second cascade species P 2 and than increases again from the second cascade species P 2 to the third cascade species P 3 , while there is an increase from the preceding cascade species to the next cascade species at the nucleus. This phenomenon is caused by the concentration differences from cell membrane to nucleus, which is larger for P 1 than for P 2 in the PC cascade. Note that the parameters were chosen to be a n b n ¼ 2, which means a twofold increase for the average concentration levels from one signaling cascade element to the next. Therefore, the spatial system can behave entirely different than the homogeneous system. The accumulation time at the membrane was much faster for γ = 1 μm/s than for γ = 0 and changed only slightly with cell size. However, for larger cells the accumulation time of the signal at the nucleus was almost independent of γ. Therefore, the difference of accumulation times at the membrane and nucleus increased with γ (also compare figure in S1 Appendix). In case of the MMC cascade the accumulation time at the nucleus for a cell with R cell = 12 μm almost doubled compared to a small cell with R cell = 2 μm, while for the PC cascade the increase of accumulation time with cell size was less pronounced.
For calculation of higher moments of the time scaling and the special case of a cell without nucleus we refer to [49]. An analysis for time scaling of a linear cascade in one spatial dimension with four elements including higher moments has been carried out in [50].

Quantifying the pathway sensitivity with respect to spatially heterogeneous signals
In the following we analyze signal transduction of heterogeneous external signals. For example, in cultures of mixed haploid yeast cell populations [40] as well as in microfluidic devices [51], Spatial modeling of the membrane-cytosolic interface in protein kinase signal transduction the external pheromone signal, which triggers a MAPK cascade, is not homogeneously distributed but forms gradients in the extracellular medium. The activated signaling cascade is spatially localized and triggers subsequent directed growth in S. cerevisiae [14] as well as S. pombe [15]. Furthermore, properties of protein-protein interactions and morphological changes can be tightly connected [52].
Therefore, we investigate the signal transduction in response to an external heterogeneous signal for same cell shapes as in Fig 3, which were a spherical cell, a rod shaped cell, a cell with one protrusion and a cell with two protrusions. These cell shapes occur for example in S. cerevisiae, S. pombe and during their response to stimulation with mating pheromone [53].
We tested the linear signaling cascade with a graded stimulus of the form where P sig 0 and P slope 0 are constants describing the basal signal strength and the slope of the signal, respectively. Here, we chose the origin of coordinates to be in the center of the cell and, therefore,x mid ¼ ð0; 0; 0Þ. In this way, we obtain an input signal gradient which increases linearly in x 1 -direction for P slope 0 > 0 and decreases linearly for P slope 0 < 0. The concentration at x mid 1 is given by P sig 0 . We tested the influence of asymmetries in cell shape in response to the graded stimulus Eq (40) and investigated the spatial distribution of the last signaling component of the MMC and PC cascade, which is P 3 . In Fig 5(A) and 5(C), the concentration profile of " P 3 on the cell membrane as well along a slice through the cell in response to a homogeneous signal is shown as control. Since the spherical cell is radially symmetric, no gradient was induced on the membrane. For the rod shaped cell, we observed a shallow gradient on the cell surface with higher concentration at the poles, the intracellular concentration profile exhibited two areas of low concentration, which were separated by the nucleus. This effect was more pronounced for the MMC cascade. For the PC cascade, the concentration was almost homogeneously distributed. For the asymmetric cell shapes with one and two protrusions, a gradient from the distal end (front) to the spherical part (back) was established in response to a homogeneous input signal. The protrusion, therefore, can be compared to a pocket in which higher concentrations of cytosolic signaling molecules are established. Mathematically this effect can be explained by the geometry dependence of the eigenfunctions of the Laplace operator [54], which can be employed to characterize the solution of the reaction-diffusion equations for a certain cell geometry. Therefore, these asymmetric cell shapes can already induce a gradient of signaling molecules from front to back.
In Fig 5(B) and 5(D), the responses to a signal with P slope 0 ¼ 0:03 mm À 1 , which is increasing in x 1 -direction, and a signal with P slope 0 ¼ À 0:03 mm À 1 , which is decreasing in x 1 -direction, were simulated and opposed to the response to a spatially homogeneous signal with P slope 0 ¼ 0. To measure the response, we define the gradient of the n-th cascade element naturally as the difference of concentrations at two points over the euclidean distance of these two points. In the case of the kinase concentrations, the gradient was computed from ðP n ðx front Þ À P n ðx back ÞÞ=jx front Àx back j. Here,x front andx back are the extreme points in x 1 -direction on the cell membrane or nucleus. Both motifs, the MMC and PC cascade, behave differently in the transduction of signal gradients. The gradient of the third cascade level P 3 along the cell membrane and the nucleus decreased for the MMC cascade with cell size for all shapes. For the PC, the gradient increased with cell size to a maximum value and then decreased for larger cell sizes, which suggests an optimal cell size for gradient detection and transmission. This effect was expected, since for small cells the concentration was almost P 3 on the cell membrane as well as along a slice through the cell in response to a homogeneous input signal P 0 100 nMμm. The spherical cell has a size of R cell = 2.5 μm with a nucleus of radius R nuc = 1 μm. All cell shapes have the same cell volume and contain a spherical nucleus of the same size. (B) Simulations for varying cell size measured as diameter d x in x-direction and three different signal slopes P slope 0 ¼ À 0:03 mm À 1 (green), 0 (orange), 0.03 μm −1 (red) were performed. The gradient grad " P 3 ≔ ð " P 3 ðx front Þ À " P 3 ðx back ÞÞ=jx front Àx back j was plotted. Here,x front andx back are the extreme points x-direction on the cell membrane or nucleus. In (C) and (D), the PC cascades was simulated for the same setup as in (A),(B). The parameters used were α 1 = α 2 = α 3 = 1 s −1 , β 1 = β 2 = β 3 = 1 s −1 , γ = 0.5 μm s −1 , D mem = 0.03 μm 2 s −1 and D cyt = 3.0 μm 2 s −1 . https://doi.org/10.1371/journal.pcbi.1006075.g005 Spatial modeling of the membrane-cytosolic interface in protein kinase signal transduction homogeneous in the cytosol and concentration differences were balanced by diffusion. However, with increasing cell size the average concentration level decreased in the cell and at the nucleus. As a consequence, also the absolute gradient decreased.
The rod shaped cell showed a stronger response than the spherical cell shape, since concentrations were higher at the poles and the cell was aligned along the gradient. Furthermore, the compartmentalization induced by the nucleus in the thin rod shaped cell had a pronounced effect on the P 3 gradient, since diffusion in the cytosol from front to back was hindered. For the cells with one and two protrusions the gradient of P 3 was strongly biased with an increase in the direction of the protrusions. Note that both motifs behave differently for the transmission of the gradient to the nucleus. While for the MMC cascade, the shape dependence was more pronounced and the gradient in the cell interior was almost decoupled from the gradient on the membrane for the asymmetric cell shapes, the PC cascade transmitted the gradient more reliably into the cell interior and the nucleus.
In summary, we observed a strong influence of cell size on localization and establishment of gradients by signaling cascade elements. For the cell with a protrusion the concentration of P 3 was higher in the protrusion than in the opposite distal end, which is the spherical part of the cell. This effect emerged due to a higher local surface to volume ratio in the protrusion region. Therefore, a larger portion of cytosolic signaling molecules, which diffuse freely in the cytosol, is phosphorylated in the protrusion part leading to a gradient from the protrusion tip to the opposite distal end of the cell. The influence of cellular asymmetries has also been investigated in [23] for gradients of the small Rho-GTPase Cdc42 during cell polarization. However, this system reacts in the opposite way, since the flux of molecules during the establishment of a polarity site is directed from the cytosol onto the membrane and, therefore, a gradient from the distal end to the protrusion is established.
These effects occur due to the different architectures of both systems. In the PC and MMC signaling cascades, we have signal transduction from the membrane to the nucleus and, therefore, a diffusive flux of activated signaling molecules from the membrane into the cytosol, while in the polarization system the flux of signaling molecules during the establishment of a polarity site is directed from the cytosol onto the membrane, which is the opposite direction. Therefore, both system respond differently to cellular asymmetries with respect to gradient formation. This interplay of both systems is especially interesting, since in many organisms a polarization system is interacting with a MAPK cascade [55,56] and might, therefore, precisely control cell shape and size.
For spherical cell shapes we furthermore investigated more complex external signal gradients, meaning heterogeneities with multiple maxima and minima (see S1 Appendix). As in [18,57], a heterogeneous signal on a sphere can be decomposed using spherical harmonics P 0 ðy; 0; tÞ ¼ In this decomposition the amplitudes of higher order, where the order is denoted by l, are generally more strongly damped than gradients or spatial heterogeneities of lower order. In this manner, the results shown here can be extended to complex spatial signals on the cell surface. We provide full analytical solutions for the MMC and PC for a sphere with excluding nucleus (see S1 Appendix).

Systems with feedback
In this section, we analyze the influence of cell size on signal transduction for an oscillating cascade consisting of two membrane-bound and one cytosolic member (MMC) and a cascade of three cytosolic elements (PC), meaning for M = 2 and M = 0, respectively. The case of a negative feedback and a constant homogeneous signal is investigated in the following. Negative feedbacks are a frequent regulation element in signaling cascades and can be induced by the dephosphorylation of upstream components by the MAPK or phosphatases [39,[58][59][60][61]. Examples are given by Tyr phosphatases, which can induce a negative feedback [47] and dual specificity phosphatases (DUSPs) [59]. Some negative feedbacks, as for instance induced in the Src-Tyr cycle, lead to oscillations on the time scale of seconds [47], while others act on much longer time scales. For instance, during the yeast pheromone response the MAPK Fus3 undergoes sustained oscillations in the range of 2-3 hours, which control the periodic formation of mating projections. In this process Sst2 acts as a negative regulator of the G-Protein signaling at the membrane, while deactivation in the cytosol is mediated by the MAPK phosphatase Msg5 [39]. A classical and most simple example of an oscillator with negative feedback and non-linear reaction terms is the Goodwin oscillatory system [62,63].
We adapt the mentioned, modified system and formulate the problem using partial differential equations by adding a diffusion term and formulating the boundary conditions accordingly to the models mentioned before. The phosphorylation and dephosphorylation rates for both models read as according to Eqs (1)-(5), respectively. The activation rate v a 1 contains the negative feedback, since a high concentration of P 3 leads to a lower activation of P 1 . We assume a constant external signal P 0 ðx; tÞ ¼ 100 nMmm.
For the first model (MMC), the deactivation with rate v d 3 takes place in the cytosol, whereas the activation occurs on the membrane and is therefore modeled as a boundary condition with v i = 0 according to Eqs (2) and (3), as P 3 is a solely cytosolic species. We assume zero-flux conditions for P 3 on the nucleus, meaning that we set the nuclear-import reaction rate = 0 (compare Eq (8)).
For the second model (PC), all species are solely cytosolic, hence the activation rate v a 1 for P 1 is a boundary condition describing the activation of P 1 on the membrane. We assume a zero-flux condition for P 3 on the nucleus ( = 0) and for P 1 and P 2 on both nucleus and membrane, meaning the whole boundary.
Both models contain non-linear kinetics as well as negative feedbacks, resulting in oscillations. Furthermore, in both models the activation rate for P 1 depends on a parameter p > 0. It is shown, e.g. in [64], that for p > 8 the ODE system consisting of three species destabilizes, and that for longer cascades, i.e. for larger N, the system becomes instable for even lower values of p > 1.
Since an analytical solution is unknown for both models, numerical methods have to be employed to solve the systems. For simplicity reasons and due to the high computational overhead we solved the systems in two dimensions, using a disc to model the cell. We used a fixedpoint scheme to solve the non-linear equations.
We chose the parameters β 1 = β 2 = β 3 = 0.125 s −1 , D cyt = 1 μm 2 s −1 , D mem = 0.01 μm 2 s −1 , R nuc = 1 μm and R cell = 2 μm. The initial conditions were P 1 = P 2 = 10 nMμm and P 3 = 10 nM for the MMC cascade and P 1 = P 2 = P 3 = 10 nM for the PC cascade. For K m = 100 nM and the feedback strength p = 10, both models oscillated as expected, Fig 6A. Therefore, in the case of a relatively small cell size of R cell = 2 μm, both spatial models behave similarly to the original model, which was formulated as a system of ordinary differential equations. An analysis of the oscillation frequencies and the mean concentration can be seen in Fig 6B. Based on previous experiments and plots, the frequency analysis was conducted after t = 200 s, when the frequency and corresponding amplitudes of the mean concentrations for both models can be assumed to be constant. Both models show a very similar behavior, Spatial modeling of the membrane-cytosolic interface in protein kinase signal transduction whereas the frequency and mean concentration are higher for the pure cytosolic model, as can be explained by the fact that the reactions do not only occur at the membrane, but everywhere in the cytosol.
In a next step we varied the cell size 1.5 μm R cell 4.0 μm and again conducted an analysis of the frequencies and corresponding amplitudes for both models. The results of the analysis for the third component of both models are plotted in Fig 6C. As pointed out before, the frequency for both models is very similar, but the amplitude of the signal is higher for the cytosolic model.
Oscillations in the first, mixed cascade model only occur for a cell size up to R cell 2.5 μm, and in the second, pure cytosolic model for a cell size up to R cell 3.0 μm. The inital oscillations die down fast and both models converge against a steady state if the cell size is chosen bigger.

Discussion
Stimulated by the progress in cell imaging and the increasing need to understand intracellular dynamics, we investigated and discussed a general approach of modeling cellular signal transduction in time and space. Signaling cascades of covalent protein modifications, such as mitogen-activated protein-kinase (MAPK) cascades and small GTPase cascades, occur in a plethora of variations [1,13,65]. The first signal component can be activated at the cell membrane by a membrane-bound enzyme such as a kinase or a guanine nucleotide exchange factor in GTPase signaling, while deactivation can occur at the membrane or in the cytosol, for instance, mediated by a phosphatase or GTPase activating protein [66]. Therefore, activation and deactivation can be spatially separated, which creates a number of different spatial arrangements and combinations in signal transduction.
We investigated signaling cascades with different spatial arrangements of signaling components. We showed that modeling of the membrane-cytosolic interface is crucial as well as the ratio of membrane area and cytosolic volume, which are both spatial properties. The results imply strong cell size and shape dependence of signal transduction within cells, which are likely to contribute to single cell variation in response to extracellular stimuli. We suggest that cells measure the cell membrane to cell volume ratio to coordinate growth and differentiation. For asymmetric cell shapes also local changes in cell volume to cell membrane ratio becomes important for intracellular signaling. Widely used time-dependent models of ordinary differential equations can naturally be extended into space by using bulk-surface differential equations. Applying this extension to a class of linear signal transduction models, we compared the assumption of a well mixed cell with two different spatial signal transduction motifs. We derived and discussed criteria that can be used to test the well-mixed assumption and showed that kinetics that connect membrane-bound species with cytosolic species naturally cause size dependence. The results are, therefore, of general importance for kinetic models of signal transduction.
Our findings have relevant biological implications. Since the signals transduced by linear signaling cascades from the cell membrane to the nucleus decrease exponentially on a length scale of a few microns, our theoretical findings suggest a strong cell size dependence in response to extracellular stimuli. Furthermore, the global cell volume to cell membrane area is important for average concentration levels. Mating projections as they occur in yeast act as pockets for signaling molecules, which can support biochemical feedbacks. Adaptations as lamellipodia in keratocytes or invaginations such as T-tubuli in myocytes can locally increase the accumulation of signaling molecules. These cellular structures are able to directly provide a feedback on signaling.
We suggest the normalized variance as a measure to quantify concentration differences and localization of signaling molecules, which can be obtained from spatially resolved microscopy data additionally to mean intensity levels of a fluorescence marker. For example, it would be enlightening to measure average concentration and normalized variance together with cell size and morphology. Interesting studies of the response in cell populations often lack the response behavior attributed to cell size and morphology. Examples range from the switch-like behavior in populations of oocytes [67] to the pheromone response in yeast cells [68,69]. Therefore, single cell data where the cell size is assigned to these measurements is needed for a faithful quantitative investigation of the pathway, to disentangle biochemical properties of protein-protein interactions and morphological properties such as size and shape of whole cells. Targeting signaling proteins by lipidation modifications such as palmitoylation, prenylation or myristoylation [2,3,5] could change the sequestration of a signaling cascade from a pure cytosolic (PC) cascade to a mixed-membrane bound (MMC) cascade. In the case of the mixed-membrane MMC the geometry and size dependence is more pronounced, since the first signaling elements are tethered to the membrane. In contrast, for the investigated PC cascade, localization and strong intracellular gradients are reduced, but depending on the kinetic parameters, the geometric information can also be better transmitted through the whole cell.
In non-linear signaling systems, the differences that we observed in the linear signaling cascade models are likely to be amplified. Non-linear kinetics can amplify gradient formation, which leads to even stronger intracellular concentration differences [70]. This also holds for absolute concentration levels that can behave in a switch-like manner depending on the kinetics [67,71]. Furthermore, higher order kinetics can amplify the accumulation time differences in different cellular locations [72], which can lead to spatial oscillations and phosphoprotein waves.
The analysis of the signaling cascade model can be extended to more complex spatial heterogeneities for example by using the Laplace series as suggested in [18,57]. With this approach localized signals arising from membrane structures like lipid rafts, septins or colocalization due to protein-protein interactions can be represented. Since these are often precursors for cell shape and organelle structures, the interplay with cell shape and morphology needs to be addressed by future research. The intrinsic geometry dependence of signaling systems has recently been shown for ellipsoidal cell shapes in the MinE-MinD system [24,73,74], but also in the yeast system [23,[75][76][77]. Recent developments of mathematical methods such as the finite element method for bulk-surface equations [19,20] as well as stability analysis techniques of these systems [23,[78][79][80][81][82] are expected to provide further insight in the behavior of cellular signal transduction.

Methods
We used the finite-element software FEniCS [83,84] to solve the arising partial differential equations in the Python programming language. The meshes were generated using the computational geometry algorithms library (CGAL) [85]. The non-linear equations were solved using a fixed-point scheme.
Supporting information S1 Appendix. The appendix contains all derivations of the analytical solutions for the steady state of the MMC and PC cascades as well as a general cascade with an arbitrary number of elements. The correspondence of the homogeneous ordinary differential equation system to the spatial MMC and PC system is established. Furthermore, an analytical expression for accumulation time of the MMC cascade is derived. (PDF)