Geometric principles of second messenger dynamics in dendritic spines

Dendritic spines are small, bulbous protrusions along dendrites in neurons and play a critical role in synaptic transmission. Dendritic spines come in a variety of shapes that depend on their developmental state. Additionally, roughly 14–19% of mature spines have a specialized endoplasmic reticulum called the spine apparatus. How does the shape of a postsynaptic spine and its internal organization affect the spatio-temporal dynamics of short timescale signaling? Answers to this question are central to our understanding the initiation of synaptic transmission, learning, and memory formation. In this work, we investigated the effect of spine and spine apparatus size and shape on the spatio-temporal dynamics of second messengers using mathematical modeling using reaction-diffusion equations in idealized geometries (ellipsoids, spheres, and mushroom-shaped). Our analyses and simulations showed that in the short timescale, spine size and shape coupled with the spine apparatus geometries govern the spatiotemporal dynamics of second messengers. We show that the curvature of the geometries gives rise to pseudo-harmonic functions, which predict the locations of maximum and minimum concentrations along the spine head. Furthermore, we showed that the lifetime of the concentration gradient can be fine-tuned by localization of fluxes on the spine head and varying the relative curvatures and distances between the spine apparatus and the spine head. Thus, we have identified several key geometric determinants of how the spine head and spine apparatus may regulate the short timescale chemical dynamics of small molecules that control synaptic plasticity.

organelle called the spine apparatus (SA), which is a protrusion of the smooth endoplasmic reticulum (ER). The spine apparatus may also be of various shapes that can change in response to signaling stimuli over days 6,[38][39][40][41] . It is generally accepted that changes in dendritic spine calcium levels, as well as localized protein synthesis, play a central role in structural plasticity, and these two processes may be influenced by the presence and shape of a spine apparatus. Experimental observations have shown that mice lacking spine apparatus (synaptopodin-deficient mice) show deficits in LTP and impaired spatial learning, thus supporting the hypothesis that the spine apparatus is involved in synaptic plasticity [42][43][44][45][46][47][48] . Furthermore, it has been hypothesized that computational modeling of the processes underlying structural plasticity could help to identify the regulatory feedback that governs the switch between LTD and LTP in ER-containing spines 49 .
Despite the emerging importance of the role of spine shape and internal organization in synaptic plasticity, the precise nature of how the physical aspects of a dendritic spine affect signaling dynamics of second messengers such as IP 3 , Ca 2+ , and cAMP remains poorly understood. In this paper, we conducted a systems biophysics study with the goal of identifying some of the design principles associated with the regulation of second messenger dynamics in dendritic spines. Specifically, using a combination of theory and computation, we sought to answer the following questions: (a) How are the spatio-temporal dynamics of second messengers within the spine affected by spine geometry? (b) How do the presence, size, and shape of the spine apparatus affect these spatio-temporal dynamics? And (c) how does different localization of the postsynaptic density (PSD, a protein-dense region in the postsynaptic membrane) affect the spatio-temporal dynamics of the second messenger? Using a minimal model to represent reaction-diffusion events within the spine, we showed that in the short timescale, the geometric features of the spine play an important role in establishing the spatio-temporal dynamics of the second messengers.

Methods
Governing equations. We developed and analyzed a reaction-diffusion model of second messenger dynamics with time-dependent flux boundary conditions. The resulting system of equations was analytically and numerically solved in simplified geometries to identify how the dynamics of second messengers are related to the geometrical parameters (see Fig. 1b). We considered a second messenger with concentration distribution C = C(x, t), where x is the vector of the spatial coordinates and t is time. In the volume of the domain, the dynamics of C are then given by the following partial differential equation (PDE): 2 where D is the diffusion constant of the species C, ▽ 2 represents the Laplacian operator in three dimensions, and τ is a time constant. τ represents a decay time constant associated with C. In the case of Ca 2+ , τ can be interpreted as the effective binding rate of rapid buffers. For IP 3 , τ represents the rate of degradation and for cAMP it represents the activity of phosphodiesterase. Our main goal in this study was to explore the solution to Eq. (1) for different geometries. Therefore, we chose a constant value of τ = 50 ms. Specialization of this model to Ca 2+ and cAMP can be found in Bell et al. 50 and Ohadi et al. 51,52 , respectively.
Boundary conditions. To completely define the dynamics of C in the domain, we need to prescribe boundary conditions on both boundaries of the domain. Since the dynamics of the receptor-mediated events at the PM and on the organelle membranes are time-dependent fluxes due to signaling reactions 53-56 , we prescribed the following boundary conditions The fluxes J PM and J ER usually depend on many nonlinear reaction terms. For the purposes of our analyses, we used biexponential functions to describe the PM flux and included a slight delay in the ER flux to write Here, the amplitude parameter γ, the time constants α and β, the amplitude of the PM-ER flux ratio ζ, and ER delay t ER are free parameters that can be fit to experimental 57 or simulation data 6 . Finally, to simulate the effect of the efflux through the spine neck, we included an outlet flux in a portion of the outer membrane defined as where K N is a constant with units of a velocity (μm/s).

Geometries.
We modeled the volume of the dendritic spine head using idealized geometries such that, to first approximation, they resemble the shape of a mature spine. We investigated three idealized geometries: spherical shells, oblate spheroidal shells, and idealized spheroidal mushroom-like geometries (see Fig. 1b-e). The dimensions of the spherical shell shape are denoted as outer radius R o and internal R i = ρR o . To highlight the combined effect of curvatures and size of the membranes, we considered oblate spheroidal shells and spheroidal mushroom-like geometry with outer eccentricity e o and inner eccentricity e i . As a result, the dimensions of the outer shell are major axis a o and minor axis b o = e o a o . The dimensions of the inner shell are major axis a i = ρa o and minor axis b i = e i a i . The mushroom-like geometry was constructed by removing the intersection between the oblate spheroidal shell (Fig. 1c) and another oblate spheroid centered at b o in the vertical axis of symmetry with minor axis b o /10 and major axis a o (Fig. 1d). In this study, the geometrical parameters ρ, e i , and e o were varied (in www.nature.com/scientificreports www.nature.com/scientificreports/ the range (0, 1]) to quantify the geometrical influence on the spatio-temporal dynamics of second messengers (C). For all simulations, we set the outer radius (and major axis) as R o = a o = 250 nm. computational tools. To assist in the derivation of the analytical solutions (Section 3.1), we used Mathematica 11.3 58 . Simulations for the dynamics of second messengers were performed using finite element methods available through the commercial package COMSOL Multphysics 5.3a with MATLAB2018a 59 . In particular, the coefficient -form PDE interface was used along with parametrized geometries. The solutions were produced with the parametric sweep utility. Wolfram Mathematica 11.3 and MATLAB2018a 60 were used for post-processing the solutions and generating the figures.

Model analysis and relevant parameters.
To compare the effect of different model geometries on the spatiotemporal dynamics of C, we defined the following metrics: • Gradient Metric: In order to quantify a local measure of the spatial gradient of the species C, we defined mean max where ||▽C|| mean (t) is the spatially averaged value of the gradient of C in Eq. 5 and ||▽C|| max is the maximum peak value. This quantity gives us insight on the localization of second messengers. When Ψ(t) tends to zero, there is no concentration gradient in the spine. • Lifetime of the gradient: We defined the lifetime of the gradient as the time, τ Ψ , needed to let Ψ(t) become lower than a threshold value Ψ th . In this work, we used Ψ th = 25% as the value beyond which the distribution of C is uniform. This quantity gives us insight into the variation of the lifetime of the gradient with respect to different geometr ies and boundary conditions.
th For clarity, the parameters used in the model are summarized in Table 1.

Results
Analytical solution. Eq. 1 is a homogeneous PDE with time-dependent BCs given by Eq. (2). To solve it, we used the method of Generalized eigenfunction expansion after formally transforming the problem into one with homogeneous boundary conditions; as a result, the PDE becomes nonhomogeneous 61 . We defined a function, the so-called reference concentration distribution w(x, t), such that it satisfies the given nonhomogeneous boundary condition, i.e.:^∇

PM PM ER ER
We considered a solution for C such that where the unknown function v is the solution of the new nonhomogeneous PDE with homogeneous BCs. This equation is obtained by substituting Eq. (9) into Eq. (1), resulting in, : 0 , and 0 (10) is a source term for v. For a homogeneous initial condition for C, the following initial condition holds for v and w, is the initial condition for v. To solve Eq. (10), the method of eigenfunction expansion is used, which consists of expanding the unknown function v in a series of spatial eigenfunctions Λ, resulting in n n n 1 Here, T n (t) are the generalized Fourier coefficients of the eigenfunctions Λ n (x). The eigenfunctions can be found using the associated homogeneous version (Q = 0) of Eq. 10. Using the method of separation of variables, v = T(t) Λ(x), we can write, where λ is a separation constant. The separation between time and space highlights the pseudo-harmonic nature of the distribution of second messenger (C). The temporal ordinary differential equation (ODE in Eq. (13)) has an exponential decay as the solution, which is affected by both λ and the timescale τ. The spatial PDE in Eq. (13) represents the Helmholtz wave equation, which is separable in 11 three-dimensional coordinate systems, and for each of them, there exists a specific family of solutions, known as the harmonic wave functions 62,63 . However, the presence of internal and external time-dependent boundary conditions and the transcendental nature of the kernel of the solution limit the possibility of finding explicit solutions and restrict the solution to an implicit form. In what follows, we further specialize this analysis for spherical domains with and without an inner boundary (Sections 3.1.1 and 3.1.2) and oblate spheroidal shells (Section 3.1.3).
Analytical solution in spherical shells with uniformly distributed BCs. We first considered a spherical shell with an inner boundary of radius R i , an outer radius of R o , and with uniformly distributed influx on both membranes (see Fig. 1c). We used a spherical coordinate system (x = (re r , θe θ , φe φ )) defined in terms of the Cartesian (x, y, z) coordinates as shown in Fig. 2a, where r ∈ [0, ∞), θ ∈ [0, 2 π), and φ ∈ [0, π] represent the radial, the azimuthal, and polar coordinates respectively 64 . Exploiting the spherical symmetry of both the domain and the boundary condition, we assumed that C(x, t), and thus v(x, t) and w(x, t), depend only on the radial coordinate r. That is, we assume that there is no dependence on the angular coordinates θ and φ. Thus, the Laplacian and normal gradient operators are: In this framework, the simplest way to define w(r, t) respecting the conditions in Eq. (8) is The homogeneous initial condition for C and Eq. (15), result in the initial condition for v as, The solution of the Helmholtz Equation (Eq. (13)) in spherical coordinates is the sum of spherical Bessel functions of the first kind (j) and the second kind (y) of zero order [61][62][63] , where a and b are integration constants. Due to the homogeneity of both the PDE and BCs, there exists a non-trivial solution for v if we impose the determinant condition 61 . The BCs in Eq. (10) now become The eigenvalues λ n are the roots of the characteristic polynomial After some algebraic manipulation and exploiting the properties of the Bessel functions, it is possible to rearrange p(λ) as is a transcendental function and the zeros can be found only numerically as showing in Fig. 2b.
Assuming that the eigenvalues λ n have been found, the expansion in Eq. (12) can be specialized as follows,  10), we obtain two nonhomogeneous first order ODEs for the Fourier coefficients that must respect the initial conditions in Eq. (22), Despite the complexity of these coefficients, we observed that the amplitude of the wave functions in Eqs (24) and (25) are implicitly related to the temporal dynamics of the biochemical reactions on the membranes of the spine (through Q), intricately coupled by the morphology (through λ n ), and reaction kinetics (through τ).
Analytical solution in spherical domains with anti-periodic BCs. We next analyzed the effect of localized influx and efflux at the outer membrane representative of fluxes at PSD and at the neck, respectively. In this context, we investigated the analytical solution of Eq. (1) in a spherical domain with no inner boundary and anti-periodic BCs at the outer domain (see Fig. 2c) written aŝ Here, J N reflects efflux through the neck. We assumed that the amplitude of the influx and efflux are the same to simplify our analysis and obtain a semi-analytical solution by exploiting the axial symmetry of the domain and BCs. As a result, C(x, t) and therefore v(x, t), and w(x, t) depend on r and φ but not θ. The Laplacian and normal gradient operators in spherical coordinates (Fig. 2a) areφ In this framework, the simplest way to define w(r, φ, t) respecting Eq. (8) and Eq. (26) is is a rectangle function. Given the homogeneous initial condition for C and Eq. (28), the initial condition for v is Using separation of variables, Λ = R(r)Φ(φ), the Helmholtz equation in Eq. (13) leads to the following two equations To find the λ mn eigenvalues we need to enforce homogeneous BCs (∇ ⋅ | = v n 0 PM ) that correspond to finding the roots of derivative of the Bessel functions such that Assuming that the eigenvalues λ mn can be found, the expansion in Eq. (12) can be written as Exploiting the orthogonality of the eigenfunctions the initial values of the generalized Fourier coefficients are given by  10), we obtain a nonhomogeneous first order ODE for the Fourier coefficients that must respect the initial conditions in Eq. (36), The solution of the above equation leads to the following expression for a mn (t), It is worth noticing that, with respect to the uniformly distributed BCs, the asymmetry of the localized fluxes introduced new angular harmonics L m in addition to the radial wave functions J m . Thus, in addition to the complex coupling between morphology of the spine and degradation kinetics, we found that the spatial distribution of the biochemical reactions on the membrane of the spines (through φ 0 ) also regulates the dynamics of second messengers.
Analytical solution in oblate spheroidal shells with uniformly distributed BCs. We next obtained the analytical solutions for the case of the confocal oblate spheroidal shell with uniformly distributed BCs (Eq. (2)). The dimensions of the outer shell are major axis a o and minor axis b o = e o a o . The dimensions of the inner shell are major axis a i = ρa o and minor axis b i = e i a i . We used the oblate spheroidal coordinates system x = (ξe ξ , ηe η , φe φ ) related to the Cartesian coordinates as shown in Fig. 2d where ξ ∈ [0, ∞), η ∈ [−π/2, π/2], and φ ∈ [0, 2π). Surfaces with constant ξ, η, and φ are oblate spheroids, confocal one-sheeted hyperboloids of revolution, and half-planes through the z-axis, respectively (Fig. 2d). To express the Laplacian in oblate spheroidal coordinates, an alternative set of spheroidal coordinates, ε 1 = sinh ξ, ε 2 = sin η, and ε 3 = φ, can be used 62-68 as follows The normal gradient to the boundary of oblate ellipsoids is given by We assume that it is possible to establish a reference function w such that the BCs (Eq. (8)) in oblate spheroidal coordinates Eq. (40) are satisfied. In this framework, the Helmholtz equation in Eq. (13) has solution in the form Λ = w 1 (ε 1 )w 2 (ε 2 )w 3 (ε 3 ), where w 1 , w 2 , and w 3 satisfy the following spheroidal wave equations respectively 62,63,65-68 ,    Here, the eigenfunctions S, P and Q are the so-called radial wave functions, spheroidal wave functions of first kind, and spheroidal wave functions of second kind respectively, while a 1 , a 2 , b 1 , b 2 , c 1 , and c 2 are integration constants. The oblate spheroidal solution can be simplified assuming the oblate symmetry of both the geometry and the BCs and thus we write Λ = w 1 (ε 1 )w 2 (ε 2 ) and μ = 0. Also in this case, the associated eigenvalues can be found by imposing the determinant condition in an implicit form. Similar to the spherical case (Eqs (24-25) the generalized Fourier coefficients (Eq. (12)) will have an exponential time decay highlighting the pseudo-harmonic nature of the solution. As a result a deviation from the spherical shape introduces more complex spatial dependence for C. In addition to the radial variation, new angular wave functions regulate the spatiotemporal dynamics of second messengers in the φ direction.

Combined effect of spine apparatus size and diffusion coefficient on the lifetime of second messengers gradient. Thus far we have elaborated on the pseudo-harmonic nature of the solution showing
how the analytical solutions can be provided implicitly. In order to visualize these solutions, we now explore the parameter space numerically. We begin with a simulation of the spatiotemporal dynamics of C in a spherical spine head with a radius r = R o and a spherical spine apparatus with a radius R i = ρR o ( Fig. 1 and Table 1). For simplicity, we start with uniformly distributed influx (J PM (t) and J ER (t)) boundary conditions on the outer and inner boundary but with no outlet (J N (t) = 0). For illustrative purposes, we chose a set of parameters for the boundary conditions (Fig. 3a) such that we can match previously reported dynamics of calcium 6 . For a diffusion coefficient D = 100 μm 2 /s and ρ = 0.5, the profile of C at a location exactly halfway between the two membranes is shown in (red line, Fig. 3b); this temporal profile is in good agreement with previous observations (see inset in Fig. 3b). Therefore, for all the following simulations, we used the time constants (α, β, and t ER ) as shown in Fig. 3a.
The amplitude parameters (γ and ζ) have been customized to the specific geometry, to avoid abnormal peaks of concentration of second messengers. A characteristic feature of the spatio-temporal dynamics of C is the lifetime of the gradient, which is affected by both the ratio between the radii (ρ) and the diffusion coefficient (D). To decode how these two quantities affect the spatiotemporal dynamics of C, we conducted the following simulations -(i) spine apparatus size was varied by changing ρ; we used three different values of ρ (0.1, 0.5, and 0.9) to capture the extreme volume changes due to small, medium, and large spine apparatus. (ii) The diffusion constant of C was varied to capture the range of intracellular diffusion from a crowded regime to free diffusion (1, 10, 100 μm 2 /s) [69][70][71][72][73][74][75] . We found that with small www.nature.com/scientificreports www.nature.com/scientificreports/ spine apparatus (ρ = 0.1), a significant concentration gradient exists in the radial direction when the diffusion coefficient is small (D = 1 μm 2 /s) at early time points but not for larger diffusion coefficients (Fig. 4a). As the spine apparatus gets larger, there is less volume available in the spine and more surface area for the internal membrane and consequently more flux through the ER. Therefore, the variability of C is smaller even for a lower diffusion coefficient (Fig. 4a-c).
To quantify the lifetime of the gradient, we calculated the norm of the gradient of C (Eq. (5)) as a function of time. Again, when the spine apparatus is small and the diffusion coefficient is small, the gradient lasts longer. An increase in diffusion coefficient for a given spine apparatus size reduces the lifetime (Fig. 4d). As the size of the spine apparatus increases, even for small diffusion coefficients, the lifetime of the gradient decreases ( Fig. 4d-f) confirming that small spine apparatus and low diffusion will result in longer time gradients of C. On the other hand, a large apparatus, even with a small diffusion coefficient will result in rapid propagation of second messengers from the PM to the ER membrane.
Given that our results are sensitive to the value of diffusion coefficients, we surveyed the literature for estimates of diffusion coefficient of different second messengers (Table 2). Based on this survey and the fact that the spine has a highly crowed environment characterized by a viscosity 5 times higher than the other cell types 76 , we chose a value of D = 10 μm 2 /s for the remaining simulations.
Size and shape of the spine and spine apparatus affect the gradient of signaling molecules. We next investigated the effect of the shape of the spine and the spine apparatus. This is particularly relevant since dendritic spines are known to have distinct shapes and their shape is associated with function 22,30,31 . Because the shape of the spine is a result of many different geometric properties, we focused specifically on curvature by modeling the spine shapes as different ellipsoidal shells. Even though this is a mathematical idealization, the resulting solutions provide insight into how curvature variations along the ellipsoids affect the harmonic functions that govern the profile of C, for uniformly spatially distributed boundary conditions. As before, the ratio between PM-ER size is controlled by ρ and now the shape is controlled by the eccentricities of the inner and outer ellipsoids e i and e o respectively ( Fig. 1c and Table 1). We conducted a systematic variation of the magnitude of these geometrical parameters and analyzed their effect on the spatio-temporal dynamics of C (Fig. 5). The influx J PM (t) was considered distributed on all the outer boundary except for a smaller portion (0.2 μm in diameter 12 ) to include the presence of an outlet flux due to the neck J N (t) (see Fig. 1c).
For a given value of e o and e i , both set to 0.5 in this case, as ρ increases, the location of maximum concentration changes from the equator to the pole. This is because of the pseudoharmonic nature of the solutions and a complex interplay among curvatures and distance between the membranes that induces a nonuniform distribution of the surface (and thus flux) per available volume (Fig. 5a). Where the outer membrane is more curved, there is more surface per volume and this is where higher levels of concentration are attained. The increase in flux per volume www.nature.com/scientificreports www.nature.com/scientificreports/ can also be caused by changing the distance between the membranes. In fact, when this distance is too small, the curvature effects became secondary. Furthermore, when there is less volume available in the spine (increasing ρ), the overall concentration of C increases. For a constant ρ and e o , increasing e i , which results in the spine apparatus tending towards a spherical geometry (Fig. 5b), C passes from a spherical-like distribution (e i = 0.1) to a localization of the maximum on the pole (e i = 1) through an intermediate situation where the maximum is confined at the equator (e i = 0.5). These results hold even in the case of constant ρ and e i and decreasing e o (Fig. 5c). Thus, the shape and size effects of the spine and spine apparatus are a result of PM shape and ER membrane shape and the relative volume enclosed. In fact, a significant difference between maximum and minimum peaks hold for tens of milliseconds, before a well-mixed condition is reached (Fig. 5d), as showed by the asymptotic evolution of the extent of gradient Ψ (Fig. 5e).
We plotted τ Ψ (Table 1) (Fig. 5f). If τ Ψ were unaffected by geometrical parameters, we would expect to see a flat plane for each value of e i . We observed that the lifetime of the gradient depends on ρ and e o in a nonlinear manner and it is not affected by the inner eccentricity. A big ER in a spheroidal −like spine (ρ → 1 and e o → 1) represents a geometrical barrier against the diffusion of C (τ Ψ  35 ms). On the contrary, a spine with a flattened head (small e o ) with a very small or absent ER (ρ → 0) has a 4-fold shorter lifetime of the gradient (τ Ψ  8 ms). We found that the paraboloidal trend of τ Ψ from simulation is in very good agreement with the one from theory, τ Ψ  L 2 /6D. Here, the length L is defined as the semi-perimeter of the ellipses with axis a i = ρR o and b o = e o R o and calculated with the approximated (see the inset in Fig. 5g). It is worth highlighting that the simulations include coupling with the timescale of the reaction, the dynamics due to the boundary conditions, as well as their distribution at the boundary of the domain whereas the analysis is based on length scale and diffusivity alone. This results in an overall extension of the lifetime of the gradient with respect to the theoretical estimation.
The results from these simulations can be summarized as follows: first, unlike spherical shells where the spatial variation is only in the radial direction, ellipsoidal shells show a spatial variation in the z and the r directions (Fig. 5a-c); second, the spatial variation of C, particularly the location of high and low concentrations of C at a given time can switch from the equator to the pole or vice versa depending on the geometry alone (Fig. 5a-c); third, the nonlinear trend of the timescale τ Ψ can be understood in terms of classical diffusion lifetime as long as the length scales are corrected for the geometries (Fig. 5d-g). Thus, these simulations predict that a deviation from spherical shape provides spines access to a more complex phase space with respect to the properties of the gradient of C and its lifetime.

Consequences of a localized input of second messengers in dendritic spines. Thus far, we have
shown that geometry by itself, in the presence of uniformly distributed boundary conditions, produces transient localization of second messengers. However, in reality, the influx at the outer boundary is not unifor mly distributed but is mostly localized to specific regions of the spine head. In fact, as shown by the reconstruction of many dendritic spines 6,12,13,40,[77][78][79][80] , the post-synaptic density, a portion of the PM rich in ionotropic receptors is localized along different positions of the membrane. Furthermore, the spine neck acts as a sink for second messengers, allowing signal transmission toward the dendritic shaft. Therefore, we next analyzed the combined effect of spine geometries, influx BCs localized to specific portions of the boundary and presence of an outlet in the neck region (Fig. 6).
We first considered a spine where the influx is confined to the upper pole region of the head. (Fig. 6a-e). In this representative case (such as ρ = 0.9, e o = 0.8, and e i = 0.1), there is a significant localization of high concentration of C close to the pole lasting for more than 20 ms (Fig. 6b-d). The gradient extends from the top to the bottom of the spine head with the spine apparatus acting as a geometric barrier for diffusion of C. The lifetime of the gradient depends on the geometrical parameters in a similar way as the uniform influx BCs, as shown by the surfaces τ Ψ (ρ, e i , e o ) in Fig. 6e.
Furthermore, when the PM influx is localized to the side of the spine (Fig. 6f-j), a big ER in a spheroidal−like spine (ρ → 1 and e o → 1) still represents a geometrical barrier against the diffusion of C. On the other hand, the lifetime of the gradient (Fig. 6j) appears to be much shorter (τ Ψ = [8~15] ms) because the influx is localized closer to the neck, providing a shorter path to the efflux boundary. From these simulations, we conclude that a localized influx on the side of the spine reduces the dependence of τ Ψ to the geometrical parameters ρ and e o when compared with the case of influx localized on the pole of the spine or that of uniformly distributed influx (see Figs 5f and 6e,j).  www.nature.com/scientificreports www.nature.com/scientificreports/ Mushroom morphology, localization of influx and variation of the size of the spine. Until now we have considered spheroidal geometry with a fixed spine size. Depending on the brain region, spine type, and species the dendritic spine dimensions vary (0.3~1.5 μm in diameter) 12,13 . Furthermore, Bourne and Harris observed that heads with diameters bigger than 0.6 μm have a mushroom morphology 37 . Therefore, we investigated how the size and shape of the spine head and of the ER (i.e. varying R o , e o , ρ, and e i ), affect the spatio-temporal dynamic of second messengers in idealized mushroom-like geometry with localized influx and efflux in the pole and neck region, respectively (see Fig. 7a).
We first maintained a fixed shape (Fig. 7a: ρ = 0.9, e o = 0.8 and e i = 0.1) and varied the spine dimension (R o = 250, 400, and 600 nm). The increase in size slightly affects the spatial distribution and the time course of C max , C min , and C mean (see Fig. 7 and insets). This because the PSD areas scale with the spine size 12 . Furthermore, the higher available volume reduced the peak value of the gradient but simultaneously elongated its lifetime (Fig. 7e). In fact, the surfaces of τ Ψ as a function of ρ, e i , and e o for the three spine dimensions showed that the lifetime of the gradient increased up to 45 ms (Fig. 7f-h). Furthermore, from the simulations we noticed that the Finally, comparing between the mushroom-like and spheroidal geometries, both with a localized influx in the pole of the spine (Figs 6 and 7) we noticed that spatial distribution and the τ Ψ surfaces are in very good agreement. Thus, all the geometrical principles discussed for the oblate spheroidal geometries hold for the idealized mushroom-like case.

Nonlinear reaction kinetics affects the spatiotemporal dynamics of C in the long timescale.
Thus far, we have only considered linear kinetics of second messengers decay within the domain that could be related to the binding of second messengers to buffers. However, signal transduction events and transport via membranes are often nonlinear, particularly in the dendritic spine 33,[81][82][83] . We extended our model to investigate whether nonlinear kinetics would change the effect of geometry and BCs. To do this, we modified the reaction-diffusion equation in Eq. (1) to include Michaelis-Menten 84 kinetics as follows where the maximum rate is V max = C max /τ and K M is the Michaelis-Menten constant. We performed the simulations with C max = 4 μM, τ = 50 ms, and K M = 2 μM, chosen for sake of representation and compared the results obtained with linear reaction kinetics (Eq. (1)). The comparison of the case previously presented (Fig. 8a-c) showed that the complexity introduced with a nonlinear reaction does not affect the spatiotemporal dynamics of second messengers in the short timescale (t < 10 ms), but affects the long-term dynamics showing a decay with different slopes. Furthermore, the localization and value of the maximum concentration do not depend on the kinetics, but rather are governed by the geometry and BCs. The kinetics affects the dynamics when a well-mixed distribution is already reached in all the cases, exhibiting different decay rates.

Discussion
Recent experimental observations have presented detailed, high-resolution images of the architecture of dendritic spines highlighting a complex internal organization 40,[77][78][79][80] . Such observations serve to highlight the role of geometry and spatial features in cellular phenomena. In this work, we used a general framework to study the effect of spine geometry including the internal organization in an idealized mathematical model with the goal of identifying some governing principles that regulate the spatio-temporal dynamics of second messengers. To do this, we developed and analyzed a general mathematical model, in which, a reaction-diffusion partial differential equation (PDE) with time-dependent mixed boundary conditions (BCs) we re analytically and numerically solved.
We arrive at the following conclusions form this work. First, the lifetime of second messengers gradients in dendritic spines depends on the intrinsic coupling between geometry of the spines and boundary conditions. These boundary conditions reflect the signaling events that take place at the membrane. Numerical simulations demonstrated how the shape of the spine governs the transient localization of peak concentration of second messengers. The lifetime of this localization depends nonlinearly on the geometry. Furthermore, we also showed that localized BCs, confined to a portion of the boundary, reduce the effect of the presence of a big and spheroidal-like ER if the influx is on the side of the spine head. Analytical investigation (Section 3.1) showed that geometry dictates the specific kind of harmonic function for the spatial distribution of second messengers. The temporal dynamics in the long time is governed by the kinetics of signaling reaction in the cytoplasm. However, this separation of temporal and spatial effects is not straightforward. The time-dependent BCs (Eq. (2)) represent both the kinetics of the membrane reactions and the curvature of the boundary (in the ∇ ⋅Ĉ n term). Therefore, time-dependent BCs represent the coupling between shape and kinetics (Fig. 9).

Increasing Head Dimension
Second, localization of the fluxes plays an important role in governing the spatiotemporal dynamics of second messengers. The lifetime of the gradient is affected by the pole vs. side localization of the spine head suggesting that the localization of the PSD plays a crucial role in spine signaling. Third, the organelle membrane plays two roles. One is to act as a diffusion barrier and the second is to act as source or sink of BCs. An emerging idea in shape regulation of signaling is the role played by organelles such as the endoplasmic reticulum and nucleus 53,85 . In the case of spines, the spine apparatus is thought to play a critical role in governing synaptic plasticity [42][43][44][45][46][47][48] . We find that the relative organization of the two membranes, PM and the ER membrane, affect the geometric landscape through both shape and boundary condition effects. We also find that the organelles can act as a physical barrier to diffusion extending the lifetime of the gradient. And finally, at short timescales, the nature of the kinetics in the cytoplasm does not alter our conclusions but kinetics play an important role in the long timescale especially in coupled cascades 3,5,86,87 . These predictions are applicable not just to spines but also to cells in general. Even though our model is simplified, it allows us to identify some common principles by which geometry can be used to alter timescales of signal transduction. The notion that shape alone can affect signal transduction (Figs 4 and 5) is a principle that is now being well-accepted in the literature 3-5,7-10 and we now extend this idea to signaling subcompartments such as dendritic spines.  . Geometric principles of second messengers dynamics: The spatio-temporal dynamics of second messengers are affected by the interplay of shape and size, biochemical signaling, and membrane reactions through boundary conditions. Our study showed that the amplitudes of the harmonic functions found in the analytical solution and dictated by the geometry, intricately depend on the boundary condition in the shorttimescale and on the reaction kinetics on long timescale.
Despite the relatively simple model presented here, we found that the spatial dependence of second messengers at short timescales was nonintuitive and exhibited a complex dependence on geometry. One potential design principle that we have identified here is that dendritic spines may be able to ensure rapid and robust signal propagation from the head toward the neck and thus toward the dentritic shaft by combining flattened head shape with very small or absence of ER, especially if the influx is localized at the pole of the spine. On the contrary, if the goal is to localize second messengers for a longer time, such as to promote synaptic plasticity, a design that allows the growth of bigger ER in a spheroidal-shaped head would help to ensure longer lifetime of the gradient. Future efforts need to consider the dynamics of both the spine and the spine apparatus during structural plasticity to incorporate mechanochemical effects.
Based on these insights, the next steps in spine systems biology can focus on specific signal transduction pathways and use reconstructions of realistic geometries to identify how the simple mathematics presented here translate into computational biology. Additionally, experimentally advances in localized imaging of second messengers molecules will be necessary to test and validate the predictions made by computational modeling. These combined efforts will enable us to extend these simple models to biologically relevant processes.