1D-3D hybrid modeling—from multi-compartment models to full resolution models in space and time

Investigation of cellular and network dynamics in the brain by means of modeling and simulation has evolved into a highly interdisciplinary field, that uses sophisticated modeling and simulation approaches to understand distinct areas of brain function. Depending on the underlying complexity, these models vary in their level of detail, in order to cope with the attached computational cost. Hence for large network simulations, single neurons are typically reduced to time-dependent signal processors, dismissing the spatial aspect of each cell. For single cell or networks with relatively small numbers of neurons, general purpose simulators allow for space and time-dependent simulations of electrical signal processing, based on the cable equation theory. An emerging field in Computational Neuroscience encompasses a new level of detail by incorporating the full three-dimensional morphology of cells and organelles into three-dimensional, space and time-dependent, simulations. While every approach has its advantages and limitations, such as computational cost, integrated and methods-spanning simulation approaches, depending on the network size could establish new ways to investigate the brain. In this paper we present a hybrid simulation approach, that makes use of reduced 1D-models using e.g., the NEURON simulator—which couples to fully resolved models for simulating cellular and sub-cellular dynamics, including the detailed three-dimensional morphology of neurons and organelles. In order to couple 1D- and 3D-simulations, we present a geometry-, membrane potential- and intracellular concentration mapping framework, with which graph- based morphologies, e.g., in the swc- or hoc-format, are mapped to full surface and volume representations of the neuron and computational data from 1D-simulations can be used as boundary conditions for full 3D simulations and vice versa. Thus, established models and data, based on general purpose 1D-simulators, can be directly coupled to the emerging field of fully resolved, highly detailed 3D-modeling approaches. We present the developed general framework for 1D/3D hybrid modeling and apply it to investigate electrically active neurons and their intracellular spatio-temporal calcium dynamics.

1. Introduction

The level of detail, that can be achieved with experimental techniques in Neuroscience is ever increasing. Intracellular organelles, protein positioning, and messenger dynamics can be recorded in spatial and temporal sequences, (Spacek and Harris, 1997; Takahashi et al., 1999; Arellano et al., 2007; Chen et al., 2008; Popov et al., 2011; Burette et al., 2012). It is obvious that processes in neurons and the brain take place in three-dimensional space and that the three-dimensionality is employed by the biological system to regulate and differentiate highly complex signaling pathways in the brain. To address the role of time-dependent, three-dimensional signal processing and investigate the role of morphology on signaling from a computational standpoint, an inevitable step is to model critical processes in the brain in three space dimensions and in time.

State of the art tools in Computational Neuroscience reduce three-dimensional problems to one-dimensional ones. For instance, simulators for the electrical signaling in neurons, e.g., NEURON (Hines and Carnevale, 2003) or Genesis (Bower and Beeman, 1997), solve a one-dimensional numerical problem in space. While this has great advantages in many applications, foremost the computational speed of the methods that allows the simulation of large network activity, the drawback is the loss of modeling the intra- and extra-cellular space of neurons, and being able to include intracellular processes in a full three-dimensional resolution. The three-dimensional organization of neurons, e.g., the filigreed geometry of the endoplasmic reticulum, (Spacek and Harris, 1997) or the intra- and inter-cellular organization of spines, (Murase and Schuman, 1999; Arellano et al., 2007; Chen et al., 2008; Tai et al., 2008; Popov et al., 2011), demonstrates the necessity for modeling intracellular processes in a highly detailed fashion in order to capture the underlying physical concepts of cellular signaling.

The need for coupling different aspects of signal processing in neurons has been recognized in the past (Kerr et al., 2008; Wils and De Shutter, 2009; Andrews et al., 2010; Cannon et al., 2010; Oliveira et al., 2010). In addition to projects focussing on the coupling of existing simulators for parallel computing architectures (Djurfeldt et al., 2010), the influence of spatial channel distribution on the electrical properties (Cannon et al., 2010) or integrating reduced intra-cellular approximations of reaction-diffusion processes, (Resasco et al., 2012; Schutter, 2013; McDougal et al., 2013a), we focus on the topic of how the three-dimensional intracellular architecture of neurons influences intra-cellular signals and how the resulting models can be efficiently solved on different computing scales. Thus, the key factors are to incorporate an accurate morphology of the neuron including the three-dimensional intra-cellular architecture (which can include active and passive organelles), model a multi-ion system described by systems of partial and ordinary differential equations (these can include the exchange mechanisms across plasma- and organelle membranes) and allow bi-directional coupling between one-dimensional membrane potential and three-dimensional intra-cellular signaling computations. The complexity of three-dimensional, detailed simulations of multi-ion systems is efficiently handled by the simulation framework uG, where Finite Element or Finite Volume discretizations and multi-grid methods are used to solve the underlying large linear equation systems on, potentially, massively parallel systems, (Heppner et al., 2013; Vogel et al., in press). Making use of uG, we propose a method to couple classical 1D-simulators for the computation of membrane potentials, with detailed 3D-simulations within an on-line 1D/3D-hybrid simulation framework.

In this paper we demonstrate this newly developed approach using the example of intra-cellular calcium dynamics coupled to the membrane potential via voltage-gated calcium channels. While 1D simulations can be carried out with classical simulators, such as NEURON or Genesis, the simulations in 3D are based on systems of partial and ordinary differential equations (PDEs and ODEs) describing distinct intracellular processes. In order to merge the two into a 1D/3D-hybrid system, two things need to be done. In the first step one needs to generate from a 1D morphology (a compartment model geometry) a surface and volume grid for numerical simulations in 3D and in the second step ensure a bi-directional coupling where the membrane potentials of the 1D simulation are mapped onto the surface of the 3D problem as boundary conditions for the numerical problem and the computed intra-cellular ion concentrations are mapped back to compute updated membrane potentials. For this we developed automated tools to compute the necessary morphology representations of a neuron and to allow bi-directional coupling.

For the neuron surface meshing in 3D, we employ a triangulation algorithm that makes use of the coordinates and radii of anatomical reconstructions, stored in e.g., hoc- or swc-format, to compute an equivalent surface mesh of the reconstructed neuron. We used TetGen (Si, 2009) to generate, from the surface grid, an intra- and extra-cellular tetrahedral volume grid. In order to associate grid nodes from the triangular surface grid with nodes from the compartment model geometry we developed V_m2uG, a tool that uses a nearest neighbor-algorithm and kd-search algorithm to perform the mapping of membrane potential data onto the triangulated surface of the 3D-neuron.

In this paper we chose a setup consisting of a multi-compartment model of a CA1 stratum radiatum interneuron taken from Katona et al. (2011) and carried out simulations for electrical signal processing in NEURON. We then used the presented methods to couple 1D and 3D models and simulated calcium dynamics in the full three-dimensional intracellular space under various parameter settings. Based on these examples, we introduce our methods and tools for 1D/3D-hybrid modeling and show how existing models and data can be incorporated into highly detailed, three-dimensional simulations.

2. Results

The methods described here couple one-dimensional electrical models and three-dimensional models for intra-cellular signaling. We chose the CA1 interneuron from Katona et al. (2011) and its neuron morphology made available on NeuroMorpho.org (Ascoli, 2006) and on ModelDB (Migliore et al., 2003). Simulations of the membrane potential dynamics in 1D, i.e., on a compartment model level, were performed with NEURON (Hines and Carnevale, 2003; Carnevale and Hines, 2006), using a standard set-up defined in the Materials and Methods, section. Since intra-cellular processes are strongly regulated by calcium, e.g., (Milner et al., 1998; Bading, 1998; West et al., 2002; Clapham, 2007; Tai et al., 2008), we chose calcium dynamics regulated by plasma membrane-located calcium channels with a given density, thus modeling effectively a channel conductance density, and a diffusion-reaction process in the neuronal cytosol as a representative of three-dimensional, intracellular signaling in neurons. 3D simulations were carried out in uG, Bastian et al. (1997); Vogel et al. (in press). Note, that this is a representative setup which is applicable to any other 1D simulations and 3D intracellular processes. The coupling of both models here occurs on the level of calcium channels—these require the membrane potential in space and time on the plasma membrane, the local intra- and extra-cellular calcium concentrations, as well as the geometry itself. In this section we will introduce the models and the simulation set-up, methods for grid generation and membrane potential mapping, and will show simulation results for the described 1D/3D hybrid simulation approach.

2.1. The 3D Calcium Model

For this study we consider a calcium model on the continuum scale, including the following components:

1. Morphology: The morphology and thus the computational domain is defined by a standard compartment model, e.g., in the hoc-format (see Supplemental Figure S1 for an example). This morphology is then mapped to an equivalent three-dimensional computational domain.

2. Membrane potential: The membrane potential is an input parameter for the calcium channel models and is computed by the 1D simulations and provided to the calcium channel models as input data.

3. Calcium channels on the plasma membrane: Based on the models by Borg-Graham (Graham, 1999), we define N-/L- or T-type calcium channels (see Materials and Methods). Channel densities can be space-dependent, thus inhomogeneous channel distribution is possible.

4. Cytosolic calcium dynamics: In this study we consider diffusion of calcium and reaction of calcium with buffers in the cytosol. The computed calcium concentrations are mapped to the 1D model to compute calcium-dependent currents.

We can formulate the above points mathematically as an initial-value boundary problem for a diffusion-reaction model. To this end, let us denote the neuron geometry as Ω, which is a compact subset of ℝ³, such that

• Ω ⊂ ℝ³ with plasma membrane boundary Γ = ∂Ω,

• Ω := Ω ∪ ∂Ω and

• Ω ∩ ∂Ω = ∅

which defines the space-time cylinder

$\begin{array}{cc}Z:=\Omega \times ]0,T[\subset {ℝ}^3\times {ℝ}^{+},\forall (x,t)\in Z:T>t>0.& (1)\end{array}$

Ω is the problem-associated domain for the 3D model, reconstructed from the original compartmental model geometry defined, in our case, by a hoc file. Furthermore let c($\overrightarrow{x}$, t) be the calcium concentration function. Then the diffusion-reaction problem can be stated as

$\begin{array}{cc}\frac{\partial c(\overrightarrow{x},t)}{\partial t}=D·\Delta c(\overrightarrow{x},t)+R(c)\text{ in }\Omega & (2)\end{array}$

$\begin{array}{cc}c(\overrightarrow{x},0)=c_0(\overrightarrow{x}),\text{ in }\Omega & (3)\end{array}$

$\begin{array}{cc}\frac{\partial c(\overrightarrow{x},t)}{\partial \overrightarrow{n}}=\Phi (V_m(\overrightarrow{x},t),\overrightarrow{x},t),\text{ on }\partial \Omega ,& (4)\end{array}$

where D denotes the diffusion coefficient for cytosolic calcium, $\Delta :={\displaystyle {\sum }_{i=1}^3{\partial }^2/\partial x_i^2}$ is the Laplace-Operator and R(c) a reaction term that depends on the calcium concentration c. Note that, the vector $\overrightarrow{n}$ denotes the direction perpendicular to the boundary surface. Equation (3) is the initial condition, i.e., a calcium distribution for t = 0 in the cytosol defined by the function c₀($\overrightarrow{x}$). Function Φ defines a Neumann flux boundary condition, regulating the calcium flux through N-/L- or T-type calcium channels, and thus also depends on the space- and time-dependent membrane potential V_m($\overrightarrow{x}$, t). Φ in this study is specified by Borg-Graham model type calcium channels, which are listed in the Materials and Methods section.

Cytosolic interaction of calcium with mobile and stationary buffers is governed by reaction equations of the type:

$\begin{array}{cc}B:=k_{on,B_{mobile}}[C{a}^{2+}]·[B_{mobile}]-k_{off,B_{mobile}}·[Ca{B}_{mobile}]& (5)\end{array}$

as well as the conservation law for the buffer concentration:

$\begin{array}{cc}[Ca{B}_{mobile}]=B_{mobile,total}-[B_{mobile}]& (6)\end{array}$

The spatio-temporal dynamics of buffering molecules are described by a diffusion-process

$\begin{array}{cc}\frac{\partial B_{mobile}(\overrightarrow{x})}{\partial t}=D_{B_{mobile}}·\Delta B_{mobile}(\overrightarrow{x})-B& (7)\end{array}$

For stationary buffers the diffusion coefficient D_Bmobile can be set to zero. The parameter values of the model equations used throughout the paper are listed in Table 1.

TABLE 1

Table 1. Parameters used for the simulations using the 1D/3D hybrid framework.

Note, that our framework is not limited to the above example of calcium/buffer dynamics. The employed multi-physics platform uG has been successfully used in a wide variety of applications, ranging from, but not limited to, simulations of skin permeability in pharmaceutical applications (Hansen et al., 2008; Nägel et al., 2008, 2009; Muha et al., 2011) to groundwater flow studies (Grillo et al., 2010) and biogas reactor modeling (Muha et al., 2012). The latter demonstrates the use of the framework for large chemical reaction systems. The 3D intra-cellular model presented here can thus be extended to include multiple ion species that are non-linearly coupled, resulting in a non-linear system of PDEs.

2.2. Components of the Hybrid Framework

As mentioned earlier, our hybrid 1D/3D framework consists of geometry matching and the coupling of 1D membrane potential simulations and 3D intra-cellular simulations. Computational grids for numerical simulations in 3D rely on a geometry-defining surface grid and a discrete representation of the computational domain, i.e., a volume grid. Approximation of unstructured domains is ideally done by a triangular surface and tetrahedral volume grid. These grids need to be generated using the compartment geometry information, provided by the 1D model as, e.g., a hoc- or swc-file. To provide an automated way for generating surface and volume grids from anatomically reconstructed geometries, we developed novel tools and combined them with existing ones, which

1. Integrate NEURON to perform 1D simulations of membrane activity in neurons.

2. Generate triangular surface grids from hoc-style neuron morphology files.

3. Map 1D simulation data of membrane potentials from compartment geometries onto the two-dimensional neuron membrane (V_m2uG) which are used as boundary conditions for membrane potential-dependent channels in the 3D model.

4. Integrate ion concentrations on the equivalent surface grid to feedback into the NEURON simulation.

The general workflow of the 1D/3D coupling is depicted in Figure 1. Using the graph-structure data contained in neuron compartment geometry files, we generate an equivalent 3D geometry with the compartment model graph as its “backbone.” This procedure is fully automated and quality controlled and will be discussed in a forth-coming paper. Since our presented framework is not limited to the grid generator used here, we would like to acknowledge the efforts of other authors addressing surface grid generation from point-diameter information, (e.g., McDougal et al., 2013b).

FIGURE 1

Figure 1. Workflow of the 1D/3D hybrid framework. The workflow can be separated into 1D and 3D modeling and simulation and into a set-up and simulation phase. In the set-up phase a general purpose simulator is used to define the one-dimensional problem (e.g., a hoc-geometry and set-up file). The 3D setup consists of generating a three-dimensional representation of the defined neuron and of specifying the cellular and intra-cellular components of the model problem. In the simulation phase 1D membrane potential simulation data is computed and mapped to the 3D framework. There, the membrane potential data is included in simulating intracellular processes, such as calcium signaling. The computed calcium data is then mapped back to the 1D problem, where it is used to compute calcium-dependent membrane fluxes.

In Figure 2 we show the in hoc-style defined neuron overlayed by the corresponding three-dimensional neuron, which can be used for the full-spatial simulation of intra-cellular processes. The quality of the volume grid influences the numerical accuracy and stability. In particular the tetrahedral angles are a means of determining the grid quality (Deuflhard and Weiser, 2011). For the grid used in this paper we measured minimal and maximal aspect ratios of the the triangular surface grid AR_tri and the tetrahedral volume grid AR_tet respectively, yielding values between 0.14 ≤ AR_tri ≤ 0.86 and 0.21 ≤ AR_tet ≤ 0.8. As a guideline one can state that aspect ratios above 0.1 result in grid elements that do not affect numerical stability (Shewchuk, 2002; Deuflhard and Weiser, 2011; Thompson et al., 2012).

FIGURE 2

Figure 2. Overlay of the multi-compartmental model defined by NEURON (red) and the equivalent tetrahedral volume mesh by uG (blue). (A) View of the whole neuron from Katona et al. (2011). Note, that close to the soma local optimization of the morphology leads to slight divergence between NEURON and uG morphology. Optimization is performed when the multi-compartment geometry model would cause unphysiological intersections of dendrites at the soma in 3D, due to the fact that all dendrites in the multi-compartment description share a common soma node. (B) Magnified views of ROIs color-coded by rectangular boxes within the overview for comparison, showing the compartment model geometry defining the backbone of the generated surface grid.

2.3. Membrane Potential Mapping—V_m2uG

In 1D compartment models the membrane potential V_m is computed in one node per cylindrical compartment. In the 3D setting each original cylinder is now represented by a segment of the triangular surface grid, thus the membrane potential from one cylinder node needs to be mapped onto a calculated cluster of nodes in the triangular surface grid.

Associating each grid point $\overrightarrow{y}$ : = (y₁, y₂, y₃) of the 3D morphology with compartment model grid points $\overrightarrow{x}$ : = (x₁, x₂, x₃) (Figure 1) and the corresponding membrane potentials V_m over the time-course of a simulation thus requires a mapping

$\begin{array}{cc}\begin{array}{l}{V}_{m}2uG:{ℝ}^3\times {ℝ}^{+}\times ℝ\to {ℝ}^3\times {ℝ}^{+}\times ℝ\\ (\overrightarrow{x},t,V_{m_x})\mapsto (\overrightarrow{y},t,V_{m_y})\end{array}& (8)\end{array}$

that assigns a membrane potential V_my to every grid point $\overrightarrow{y}$. V_my is set to the potential V_mx associated with the nearest neighbor $\overrightarrow{x}$ : = (x₁, x₂, x₃) of grid point $\overrightarrow{y}$ with respect to all points contained in the hoc-file in every timestep t of the simulation. The distance measure for calculating the nearest neighbor can be any Minkowski metric and is interchangeable. In this study we used Algorithms 1 and 2 and chose the euclidean metric, or L₂-norm, such that

$\begin{array}{cc}\text{dist}(\text{p},\text{q}):={\left({\displaystyle \sum _{i=0}^{d-1}{(p_i-q_i)}^2}\right)}^{1/2}& (9)\end{array}$

Because the geometry information provided by the hoc-file may be very large, i.e., the number of provided points n is large, we implemented and tested two strategies for determining the nearest neighbor. The first is an exact solution of the problem by space partitioning with multi-dimensional binary search trees, k-d trees with a search complexity (n · log(n)) (see Bentley, 1975 as well as Figure 3 for an example), the second strategy is the exact solution by pairwise comparison with a search complexity (n²). Note that we consider a three-dimensional example here, but the mapping process is dimension-independent and can be done for arbitrary dimensions by means of the underlying C++ library for multi-dimensional (binary) search trees (Mount and Arya, 2012), a library that is used in V_m2uG.

FIGURE 3

Figure 3. Basic operations. (A) General description of the nearest neighbor search task in 2D/3D: Let S be a set of #n points with dimension d: = 2, 3. Let q be a query point. The task is to find a point r (green) in the set S of all points (black), which is closest to the query point q [center of sphere in (A) on the right] within a given bounding box, i.e., within a circle (black) or sphere. (A) Pairwise comparison of all points in a radius search. B: Space-partitioning approach using k-d trees (multi-dimensional binary search trees). (C) An exemplary multi-dimensional binary search tree of dimension d = 2 (enhanced and modified according to Wikipedia, 2012), representing the space partition in (B).

ALGORITHM 1

Algorithm 1. Linear interpolation

ALGORITHM 2

Algorithm 2. Bilinear interpolation

2.4. Algorithms

The pairwise comparison is an exhaustive search that requires no additional data structure but (n · m) iterations which is impractical if the number n of provided points on the grid is large, and for the number of corresponding points in the hoc file holds n ≈ m, thus arriving at quadratic runtime complexity (n²). For large n, we therefore use the method of space partitioning (k-d trees, Figures 3B,C) which leads to an improved average runtime complexity for a query of $\text{log}\left(\frac{n}{d}\right)$, with d being the space dimension. With this approach we reach super-linear runtimes after the initial tree has been built. In contrast to the (n²) method, we need an additional structure, which is justified because of the overall speedup compared to the naïve approach. The overall runtime is dominated by the utilized search algorithm used during the initial build of the tree and the balancing of the tree itself, yielding a complexity (n · log(n)) (Wald and Hvran, 2006; Cormen et al., 2011). The worst case arises if n « 2^d is satisfied and the runtime complexity increases per query to (Lee, 1997). Yet, the curse of dimensionality is not an issue in the presented application, since our geometry is defined in three-dimensional space.

V_m2uG offers linear or bilinear membrane potential interpolation (Algorithms 1 and 2). This can be useful if the compartmental representation of the neuron is very coarse and mapping distances become large, i.e., 3D surface grid points lie far away from the computed nearest neighbor on the compartment geometry. Interpolation is then carried out over n ∈ ℕ⁺ pseudo nearest neighbors to compute an optimal approximated value for the membrane potential which is then assigned to the given grid point.

2.5. Workflow and Bidirectional Coupling

We use a direct data coupling mechanism termed a type A problem in Heterogenous Multiscale Modeling, see (Weinan et al., 2003), utilizing an online algorithm within uG. The NEURON library is used as a plugin/shared library within the uG project. The simulation workflow is defined in uG by a lua-script (Ierusalimschy et al., 2013) (see Figure 4 for the script used in this paper). We provide fully flexible control over NEURON within uG by means of a custom API or by directly including the hoc-script (which is then executed by the HOC language interpreter). For the latter we developed an additional C++ wrapper.

FIGURE 4

Figure 4. Script file illustrating the setup of the NEURON and uG model, here divided into a NEURON setup phase, a uG setup phase, chemical model specifications and the simulation run control. One can make use of the hoc interpreter features within the lua script files that define the uG-setup. Note that one can specify the geometry and stimulation protocol in distinct files and combine this with using the hoc interpreter from uG to refine the setup by statements of the form: HocSetup:execute(“command string”).

The workflow defined in the lua-script performs the steps illustrated in Figures 1, 4. The numerical procedures are defined individually for NEURON and uG (see Materials and Methods), which is necessary since the model equations for the 1D and 3D model are of different types. The only coupling requirement is that both simulations can synchronize their data at specific time points, which is guaranteed by our framework.

Data exchange between 1D and 3D model is bidirectional, where membrane potential data is passed from 1D to 3D via the mapping algorithm presented in the previous section. The computed ion concentrations in the 3D model are then passed from 3D to 1D in the following way:

Consider a cylindrical compartment C_i of the 1D model. By ∂(C_i ∩ ∂Ω) we denote the portion of the cell surface ∂Ω that is associated with compartment C_i. We can then compute the intra-cellular calcium concentration for C_i as

$\begin{array}{cc}\begin{array}{l}[C{a}^{2+}]C_i=\frac{\Vert \partial C_i\Vert }{\Vert \partial (C_i\cap \partial \Omega )\Vert }{\displaystyle {\int }_{\partial (C_i\cap \partial \Omega )}{u}_{C{a}^{2+}}}(\overrightarrow{x})\\ \text{ for }\overrightarrow{x}\in \partial (C_i\cap \partial \Omega )\text{ and }\forall i=0,\dots ,N\end{array}& (10)\end{array}$

where u_Ca²⁺($\overrightarrow{x}$) is the calcium concentration in node $\overrightarrow{x}$ of the 3D grid and N is the number of compartments representing the 1D geometry. ‖.‖ symbolizes the size of ∂ C_i and ∂(C_i ∩ ∂Ω) respectively.

The factor $\frac{\Vert \partial C_i\Vert }{\Vert \partial (C_i\cap \partial \Omega )\Vert }$ accounts for the fact, that the surface size of ∂(C_i ∩ ∂Ω) and the surface of C_i are not necessarily identical. Computing the calcium fluxes in uG and in NEURON using the values computed by Equation (10) shows good agreement between the two models, see Supplemental Figure S2 as well as Supplemental Figure S3.

2.6. Simulations Using the 1D/3D Hybrid Method

We applied the described method for coupling 1D and 3D simulations to a CA1 interneuron from Katona et al. (2011), taken from the database ModelDB (Migliore et al., 2003). We set up a NEURON simulation with stimulation protocol (see Materials and Methods) and generated a three-dimensional geometry for calcium simulations that use the membrane potential data from the NEURON simulation (Figure 5) and feeds back the computed concentrations (see previous section). An overlay of hoc- and ugx- (the uG grid format) geometry shows good agreement between the two.

FIGURE 5

Figure 5. Geometry specifications. (A) Three-dimensional neuron from Katona et al. (2011) used in the simulations. The box highlights the region of interest used in (B,C). (B) Triangular surface grid of a dendritic segment including an intracellular obstacle (yellow). The arrow in (A) indicates the position of the obstacle. (C) Triangular surface grid of a dendritic segment without the intracellular obstacle. The surface grid can remain unaffected by the insertion of an intra-cellular object.

The scientific gain from 3D simulations lies in a physically detailed representation of morphology and the underlying biophysical processes in neurons. These processes are regulated by space and time-dependent variations of parameters like channel densities, cytosolic and organelle architecture, diffusivity and of course cellular and organelle morphology. Previous work has demonstrated the importance of investigating the spatial organization of ion channels and their stochastic behavior, (Oliveira et al., 2010; Brandi et al., 2011), in part using a multi-scale modeling approach. Brandi et al. (2011) for instance present in their abstract a platform for combining NeuroRD (Blackwell et al., 2013) and MOOSE (Ray and Bhalla, 2008). The framework employed here can be viewed as a natural extension of previously published multi-scale approaches, with a strong focus on resolving the intra-cellular space with additional accuracy, rather than interpreting the intra-cellular space as a homogenous space for reaction-diffusion processes. Obstacles or active organelles, such as calcium-relevant stores (e.g., the endoplasmic reticulum or mitochondria), can be easily integrated into our hybrid framework. In addition to this, our framework aims at simulating entire neurons up to small networks with a high level of cellular and intra-cellular detail, rather than investigating e.g., the single molecule level by Smoldyn Andrews et al. (2010). To accomplish this, we model at the continuum scale and make use of the massively scalable code uG, (Heppner et al., 2013), demonstrated to perform ideal weak scaling up to 64 k processes on the Hermit Supercomputer. In Table 2 we show ideal weak scaling for the presented simulations on our in-house cluster (see Materials and Methods), solving a 3D problem with 256,000 grid points in under 1 min on 64 processors. In comparison to other simulators, such as VCell or MCell, the VCell website lists a 2D benchmark with a problem size of 5000–10,000 grid points and a PDE-system with 4–5 PDEs being solved in roughly 20 min. In Balls and Baden (2004), MCell scaling studies up to 64 cores show a scaling efficiency of 85–92%, solving a benchmark problem also in roughly 20 min. Since the underlying models of the compared simulators differ, a direct comparison is not necessarily valid. The presented hybrid method, based on a continuum model is applicable ideally for whole- and multi-cell simulations, while e.g., MCell uses particle-based methods that so far do not go beyond spine simulations, i.e., microdomain simulations.

TABLE 2

Table 2. Scaling.

In order to demonstrate some advantages of a 1D/3D hybrid approach, we defined a standard 1D/3D simulation protocol (see Materials and Methods) and then varied parameters, such as calcium diffusivity (Allbritton et al., 1992), voltage-gated calcium channel (VGCC) densities (Pumplin et al., 1981; Eggermann et al., 2011), and added intra-cellular obstacles to the cytosol.

Figure 6 and the supplemental movies show calcium simulations in a region of interest using the 1D/3D hybrid method, with the calcium diffusion set to 100 μm²/s, a VGCC-density of 1000 μm⁻² and no obstacle present in the cytosol (see Figure 5 for the neuron geometry used in the simulations). This example shows how the NEURON-calculated membrane potentials are mapped to the 3D model, where they regulate VGCCs and detailed calcium dynamics can be simulated in full three-dimensional space. The calcium concentrations are then mapped back to NEURON according to the previous section.

FIGURE 6

Figure 6. Time series for a simulation of intra-cellular calcium signaling with a calcium diffusion coefficient set to 100 μm²/s, VGCC density of 1000 μm⁻² and no obstacle present in the dendrite. The basal calcium concentration is set to 100 nM. (A–F) show the propagation of calcium in a dendrite, triggered by the the stimulation protocol listed in Table 1. Membrane potentials were computed in NEURON and mapped to the plasma membrane in order to compute the VGCC-dynamics regulating the calcium exchange. Calcium concentrations were then mapped back to compute the calcium-dependent membrane potential in NEURON. (A–F) covers 1s of real time. Time points shown are at 0, 0.01, 0.02, 0.03, 0.1, and 1 s. Note, the maximum single-channel conductance was set to 60 pS according to Graham (1999). (G) Highlighted measuring points close to the plasma membrane (pm) in between and at the center used in (H). (H) Evolvement of the calcium concentrations at three different points over a period of 20 ms. Note that the VGCCs close at the peak amplitude and the calcium profile quickly adjust to a uniform value based on the intra-cellular diffusion.

2.6.1. Varying the diffusion coefficient

We then varied the diffusion coefficient for intracellular calcium between 10 and 100 μm²/s, which is within the experimentally observed ranges, (Allbritton et al., 1992), while the VGCC density was set to 1000 μm⁻² (Figure 7A). Figure 7B shows a linear dependence to a global variation of the cytosolic calcium diffusion coefficient. Note, that the diffusion of calcium does not have to be isotropic, but can have different diffusion properties in the three space-dimensions. The diffusion coefficient then needs to be formulated as a diffusion tensor

$D:=\left(\begin{array}{ccc}{d}_x& 0& 0\\ 0& d_y& 0\\ 0& 0& d_z\end{array}\right),$

which can be defined in the uG-workflow (Vogel et al., in press).

FIGURE 7

Figure 7. Parameter variation. (A) Variation of the diffusion coefficient regulates the intracellular calcium profile (using a constant VGCC density). The concentration profiles computed by the hybrid framework are validated in Figure S3 using an alternative, reductionist approach. (B) Quantification of the peak amplitudes in (A), R² = 0.9461 and p = 0.003, showing a linear dependency between the diffusion coefficient and the peak amplitude. (C) Density variation of VGCCs sharpens the concentration profile in the considered dendrite (using a constant diffusion coefficient). (D) Quantification of the peak amplitudes in (C), R² = 0.99 and p = 0.018 showing a linear dependency between VGCC density and peak amplitudes. Calcium profiles are validated by calcium imaging data, Helmchen et al. (1996); Kits et al. (1997); Maravall et al. (2000). VGCC densities are varied according to documented physiological intervals up to 10000/μm², as reported in Eggermann et al. (2011). (E,F) Placing obstacles of various sizes in the dendrite, from 0 to 100 % (see Table 3) affects the locally available calcium concentration, (E) before the obstacle and (F) behind the obstacle.

2.6.2. Varying channel densities

Next, we decided to vary the VGCC densities on the plasma membrane. We placed N-type VGCCs (see Materials and Methods) homogeneously along the dendrites at densities between 200 and 1000 μm⁻² (Figure 7C). Due to an increased channel density we observe changes in the rise and decay times, as well as different peak amplitudes in the calcium profile. In Figure 7D we see that the peak amplitude of the intra-cellular calcium signal depends linearly on the VGCC density. Global and time-dependent changes in the VGCC distribution of the neuron can thus finely regulate the intra-cellular calcium code.

2.6.3. Including intra-cellular obstacles

Where in a 1D model one needs to average the intra-cellular space and approximate morphologies, the 3D approach allows a full and detailed investigation. To demonstrate how intra-cellular objects can be included into a simulation, we investigated calcium dynamics in the observed dendrite with and without an intra-cellular obstacle. For this we placed an obstacle (this could be e.g., a mitochondrion or endoplasmic reticulum) in the center of a dendrite (see Figure 5), defined as a cylinder that is 0.14 μm long and has a variable circumference C_obstacle, embedded in a dendrite of average local circumference C_neurite. We then define the hindrance factor H of the intracellular obstacle with respect to the surrounding cylinder for that particular section as:

$\begin{array}{cc}H:=\frac{C_{obstacle}}{C_{cylinder}}\in [0,1]& (11)\end{array}$

The obstacle sizes used in our simulations are listed in Table 3, which correspond to a dendritic occupancy between 0 and 100% (100% corresponds to a full blockage in the the dendrite). We chose the obstacle location to be right before the most distal bifurcation point of a certain neurite which emerges—and can be visually traced in Figure 5A—from the left most bottom of the soma itself. Note, that the obstacle is placed in the middle of the considered neurite in the present use case. Intra-cellular organelles, occupying parts of the cytosol, affect the intra-cellular resistivity. To account for this in our 1D model, one can either directly modify the axial resistance by changing the intra-cellular resistance or by changing the diameter of the cylinders that harbor an organelle. The axial resistance R_a in a cylinder with length l and diameter d is computed by:

$\begin{array}{cc}{R}_a=R_i·\frac{4l}{d^2\pi }& (12)\end{array}$

where R_i is the intra-cellular resistance. The effective diameter can be calculated for all affected compartments by:

$\begin{array}{cc}{d}_{new}=d_{old}-d_{obstacle}& (13)\end{array}$

TABLE 3

Table 3. Hindrance factors.

Our framework automatically checks which compartments are occupied by organelles and adjusts the diameters accordingly. We observed that local intracellular signaling is affected by obstacle size and position (Figures 7E,F). In more complex situations, e.g., space occupancy by mitochondria, endoplasmic reticulum or spine apparatuses in dendrites or spine necks, the inclusion of these functional organelles in a three-dimensional model might be necessary for explaining the intracellular dynamics of the neuron. Active intra-cellular organelles can be easily added within the presented framework. The user needs to define the geometry of the organelle (which can be a detailed three-dimensional reconstruction) and the biological exchange mechanisms across the organelle membrane (e.g., IP3-receptors, ryanodine receptors, SERCA pumps, sodium-calcium exchangers etc.).

3. Discussion

Modeling and simulating cellular and network processes are constrained by the complexity of the computational problem and the availability of experimental data to validate the models. Experimental techniques are evolving, shedding new light on the filigreed functioning of neurons and networks. Spatial and time resolutions are becoming fine enough to resolve positions of single proteins and receptors at synapses and the intra-cellular domain can be investigated at increasing levels of detail. On the other hand computational complexity increases when a model includes biological and spatial detail, resulting in a highly detailed three-dimensional model of neurons. Computational complexity further increases when not only one neuron, but entire networks need to be modeled.

To cope with network complexity, several approximation methods for cellular and network function have been developed over the last decades. These approximations yield zero- or one-dimensional models in space, describing neurons with point- or compartment neuron models. Numerous specialized and general purpose neuron simulators incorporate these basic abstraction techniques (e.g., Bower and Beeman, 1997; Hines and Carnevale, 2003; Gewaltig and Diesmann, 2007) and allow the user to simulate large and complex neural network structures. These models and simulators were developed and evolved not only around the limitations of computational power, but also around the availability of experimental data to validate the models. The evolution of experimental methods, computational resources and numerical mathematics motivates the development of neuron models that include greater physical detail. This brings us to three-dimensional models that include the detailed morphology of neurons, either by modeling physical processes using stochastic (e.g., Andrews et al., 2010; Oliveira et al., 2010) or deterministic models. We make use of the latter approach with the additional feature of being able to integrate obstacles and intracellular organelles, which yields models formulated on the continuum scale by means of systems of partial differential equations. Solving these equations numerically requires advanced mathematical tools. One general purpose platform is uG (Bastian et al. 1997), which we used in this paper. uG offers numerical discretization methods and efficient solvers for systems of partial differential equations on highly unstructured grids and runs on massively parallel systems (Heppner et al., 2013; Vogel et al., in press). As shown in the past, these numerical tools are applicable and efficient on a cellular level or in a networks with relatively small numbers of neurons (Xylouris et al., 2007; Wittmann et al., 2009; Xylouris, 2013).

Given the limitations in computing resources, it is currently not feasible to model and simulate large networks of neurons with full single cell, three-dimensional detail. To make use of the advantages of highly detailed three-dimensional models and to cope with the complexity that comes with modeling large networks of neurons, we developed a method to couple state of the art general purpose neuron simulators, e.g., NEURON, with three-dimensional models of single neurons. This 1D/3D hybrid method includes automated tools to either reconstruct three-dimensional neuron morphologies from raw microscopy data (Broser et al., 2004; Queisser et al., 2008; Jungblut et al., 2011), from anatomically recorded data (Wolf et al., 2013) or from graph-structure morphologies as used, e.g., in the NEURON Simulator in the form of hoc-files or swc-files.

In this paper we introduced a framework for geometry and membrane potential and intra-cellular ion concentration mapping between 1D simulations and the equivalent 3D model. For this we used the NEURON simulator to compute electrical signals in a compartment neuron and uG to simulate intracellular calcium dynamics in 3D. With this approach we were able to exploit the modeling and computational advantages that a general purpose simulator for large networks brings, as well as the necessary tools to investigate intra-cellular (or even extra-cellular) processes on a very fine scale. By including the detailed morphology of neurons—which can be subject to temporal adaption due to neuronal activity (Silver et al., 1990; Korkotian and Segal, 1999; Muller et al., 2002; Van Aelst and Cline, 2004; Hayashi and Majewska, 2005; Tada and Sheng, 2006; Wittmann et al., 2009; Kanamori, 2013)—the interplay between cellular/organelle morphology and cellular function can be systematically investigated. Using defined stimulation protocols in NEURON, we ran 1D-simulations mapping the results onto the three-dimensional morphology as boundary conditions for a cellular calcium model and vice versa. The calcium model consisted of a local distribution of N-type voltage-gated calcium channels, regulated by the membrane potential and intracellular diffusion and reaction of calcium ions with a mobile buffer. After presenting the functionality of the 1D/3D hybrid method on a reference parameter set, we varied the density of voltage-gated calcium channels, the calcium diffusion coefficient and introduced intracellular obstacles. The results show that for one, the presented method is designed in a general fashion and is thus applicable to a broad range of neurobiological questions and that the effect of intra-cellular obstacles, locally changing channel densities and cytosolic diffusivity have a substantial effect on intracellular signals.

For future research, problem-specific models can be used within the presented framework. For instance, not only single neurons, but entire networks could be simulated on the 1D level, where a small subset of “key neurons” in the network can be resolved in full three-dimensional detail (Xylouris et al., 2007). Furthermore intra-cellular models for different ionic species, detailed models of intra-cellular organelles, such as mitochondria and endoplasmic reticulum calcium exchangers, protein synthesis and synapses could be included on the 3D scale.

4. Materials and Methods

4.1. The 1D Electrical Model

To simulate the electrical activity of the model neuron, we used the classical Hodgkin-Huxley equations (Hodgkin and Huxley, 1952) including sodium, potassium, calcium and a leakage current:

$\begin{array}{cc}\begin{array}{l}{i}_m=c_m\frac{\partial V}{\partial t}+n^4{\tilde{g}}_K(V-E_K)+m^{3}h{\tilde{g}}_{Na}(V-E_{Na})\hfill \\ +{\tilde{g}}_L(V-E_L)+I_{Ca},\hfill \end{array}& (14)\end{array}$

where V_m is the membrane potential, c_m the membrane capacity, $\tilde{g}$_Na, $\tilde{g}$_K and $\tilde{g}$_L are the maximum single channel conductances for sodium and potassium channels, as well as a leakage current. E_K, E_Na and E_L are the reversal potentials. I_Ca is the calcium current specified by Equation (26). The time and voltage-dependent gating variables n, m, h determine the gating behavior of the channels and are computed by the ordinary differential equation

$\begin{array}{cc}\frac{\partial x}{\partial t}=\frac{x_{\infty }-x}{{\tau }_x}& (15)\end{array}$

where x = n, m, h. The parameters x_∞ and τ_x in the equations above are computed by

$\begin{array}{cc}{x}_{\infty }=\frac{{\alpha }_x}{({\alpha }_x+{\beta }_x)}& (16)\end{array}$

and

$\begin{array}{cc}{\tau }_x=\frac{1}{{\alpha }_x+{\beta }_x}+{\tau }_{x,0}& (17)\end{array}$

where τ_{x, 0} = 0 for x = n, m, h. The values for α_x and β_x are derived from (Hodgkin and Huxley, 1952):

$\begin{array}{cc}{\alpha }_n=\frac{0.01(V+10)}{e^{\frac{V-10}{10}}-1}& (18)\end{array}$

$\begin{array}{cc}{\beta }_n=\frac{e^{\frac{V}{80}}}{8}& (19)\end{array}$

$\begin{array}{cc}{\alpha }_m=\frac{0.1(V+25)}{e^{\frac{V+25}{10}}-1}& (20)\end{array}$

$\begin{array}{cc}{\beta }_m=4e^{\frac{V}{18}}& (21)\end{array}$

$\begin{array}{cc}{\alpha }_h=0.07e^{\frac{V}{20}}& (22)\end{array}$

$\begin{array}{cc}{\beta }_h=\frac{1}{e^{\frac{V+30}{10}}+1}& (23)\end{array}$

Solving a multi-compartment model where for each cylindrical compartment an equation of the type

$\begin{array}{cc}\frac{d{V}_m}{dt}=-i_m+I_{electrode}& (24)\end{array}$

is computed using the NEURON simulator. Numerically, the arising sets of ordinary differential equations are solved with the Crank-Nicolson method implemented in NEURON. We chose the CA1 stratum radiatum hippocampal interneuron published in Katona et al. (2011) and ran the following stimulation protocol:

We inserted a current clamp at the mid of the soma and stimulated with amplitude of 0.1 nA with timestep dt = 0.1 ms for a total time of 10 ms. We then simulated 1 s of membrane potential propagation in NEURON evoked due to that stimulation.

4.2. Voltage-Gated Calcium Channel Models

In order to include voltage-gated calcium channels, we used the Borg-Graham model described for different ion channel types in detail in Graham (1999). In our setup we included N-type calcium channels gates, though the implementation allows us to also include other types (e.g., L-/T-type gates). We define a mapping of the calcium fluxes calculated by the VGCC-model onto the 3D morphology:

$\begin{array}{cc}\text{Borg Graham}:{ℝ}^3\times {ℝ}^{+}\times ℝ\to ℝ:(x_b,t,v_{m_b})\mapsto F& (25)\end{array}$

The dynamics of calcium ionic fluxes are described in the following way (Graham, 1999) and are computed inside the 1D-simulation loop:

$\begin{array}{cc}{I}_{Ca}(V,t,{[Ca]}_i,{[Ca]}_o)=G(V,t)F(V,{[Ca]}_i,{[Ca]}_o)& (26)\end{array}$

G describes the properties of the gating functions, in particular the open and close state-probabilities of the channels and is specified in Equation (28). F can be computed using the GHK model, yielding

$\begin{array}{cc}\begin{array}{l}F(V,[Ca]i,[Ca]o)\hfill \\ =p_{Ca}\frac{z^2{F}^{2}V}{RT}\left[\frac{{[C{a}^{2+}]}_i-{[C{a}^{2+}]}_o\exp \text{}(-z\frac{FV}{RT})}{1-\exp \text{}(-z\frac{FV}{RT})}\right]\end{array}& (27)\end{array}$

where p_Ca is the permeability and for N-type calcium channels is set to the value 10⁻⁸ cm³/s. F denotes the Faraday constant, R the gas constant and T the temperature. For the VGCCs used in the simulations, the gating function G is defined as

$\begin{array}{cc}G(V,t)=k(V,t)l^2(V,t)& (28)\end{array}$

Here, k and l fulfill the ordinary differential equation (15) and eqs. (16, 17), using the following values:

$\begin{array}{cc}\begin{array}{l}{\alpha }_x(V)=K_x\exp \left(\frac{z_x{\gamma }_x(V-V_{1/2,x})F}{RT}\right)\\ {\beta }_x(V)=K_x\exp \left(\frac{-z_x(1-{\gamma }_x)(V-V_{1/2,x})F}{RT}\right)\end{array}& (29)\end{array}$

The parameters V_1/2, z, γ, K, τ₀ also depend on the particular channel and are documented in Graham (1999). We set the values accordingly:

$\begin{array}{l}{V}_{1/2,k}:=-21\text{mV}\hfill \\ z_k:=2\hfill \\ {\gamma }_k:=0\hfill \\ K_k:=1.7\text{ms}\hfill \\ {\tau }_{k,0}:=1.7\text{ms}\hfill \\ V_{1/2,l}:=-40\text{mV}\hfill \\ z_l:=1\hfill \\ {\gamma }_l:=0\hfill \\ K_l:=70\text{ms}\hfill \\ {\tau }_{l,0}:=70\text{ms}\hfill \end{array}$

It is sufficient to solve Equation (15) using the first order explicit Euler method, in which the required membrane potential is taken from the membrane potential mapping (V_m2uG) in each time step.

4.3. Numerical Approach

The calcium model used in the presented simulations is based on a cytosolic diffusion equation with boundary conditions that include the plasma membrane calcium exchange mechanisms, i.e., Neumann boundary conditions representing an ionic flux over a given surface area of the plasma membrane. This can be stated as:

$\begin{array}{cc}\nabla ·(D\nabla u(\overrightarrow{x},t))=\frac{\partial }{\partial t}u(\overrightarrow{x},t)\text{ in }\Omega \subset {ℝ}^3& (30)\end{array}$

with boundary conditions

$\begin{array}{cc}\frac{\partial u}{\partial n}=g(Channels)& (31)\end{array}$

where u denotes the intracellular calcium concentration, D the diffusion tensor (a 3 × 3-Matrix) and g(Channels) the calcium flux depending on the channels defined in the simulation (see previous section). To solve Equation 30, we discretized in space by means of first order finite volumes, treating the time derivative by an implicit Euler method.

For this we discretized the plasma membrane geometry with a triangular surface grid and the volume of the cell with a tetrahedral volume grid. In the following we denote the discrete volume of the neuron as Ω and each arising control volume box as B_i, i = 1, …, n. Therefore Ω = ∪ⁿ₁B_i. The Finite Volume method (cf. Eymard et al., 2000) defines an ansatz space V_h ⊂ H¹(Ω) which is spanned by piecewise linear ansatz functions Λ_i($\overrightarrow{x}$_i) = δ_ij and is used as the approximation for the solution u. The approximation of u is then defined as the linear combination of these ansatz functions, resulting in a system of linear equations

$\begin{array}{cc}{\displaystyle \sum _j^n{A}_{ij}}{u}_j=b_i& (32)\end{array}$

where A is the stiffness matrix. To solve the linear system of equations, we applied a geometric multigrid solver (smoother ILU) as a preconditioner along with BiCGstab as the base solver for the linear part (inner loop) and a Newton solver for the non-linear part (outer loop), see e.g., (Hackbusch, 1985, 2003). The numerical setup and simulations were carried out in uG [see Bastian et al. (1997); Vogel et al. (in press)].

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

Numerous gratitude is given to Andreas Vogel pointing out important properties of the solvers and their implementation in uG, as well as Anton Vadimovich Chizhov of the Physical-technical Institute Ioffe (Russian Academy of Sciences) in St. Petersburg, Russia for providing data for validation of the Borg Graham Channel model (Graham, 1999) used in conjunction with our hybrid method.

Supplementary Material

The Supplementary Material for this article can be found online at: http://www.frontiersin.org/journal/10.3389/fninf.2014.00068/abstract

Supplemental Movie 1. Simulation of 200 ms real time, perspective 1.

Supplemental Movie 2. Simulation of 200 ms real time, perspective 2.

References

Allbritton, N., Meyer, T., and Stryer, L. (1992). Range of messenger action of calcium ion and inositol 1,4,5-triphosphate. Science 258, 1812–1815. doi: 10.1126/science.1465619

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Andrews, S. S., Addy, N. J., Brent, R., and Arkin, A. P. (2010). Detailed simulations of cell biology with Smoldyn 2.1. PLoS Comput. Biol. 6:e1000705. doi: 10.1371/journal.pcbi.1000705

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Anwar, H., Hepburn, I., Nedelescu, H., Chen, W., and De Schutter, E. (2013). Stochastic calcium mechanisms cause dendritic calcium spike variability. J. Neurosci. 33, 15848–15867. doi: 10.1523/JNEUROSCI.1722-13.2013

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Arellano, J. I., Benavides-Piccione, R., DeFelipe, J., and Yuste, R. (2007). Ultrastructure of dendritic spines: correlation between synaptic and spine morphologies. Front. Neurosci. 1, 131–143. doi: 10.3389/neuro.01.1.1.010.2007

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Ascoli, G. A. (2006). Mobilizing the base of neuroscience data: the case of neuronal morphologies. Nat. Rev. Neurosci. 7, 318–324. doi: 10.1038/nrn1885

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Bading, H. (1998). Transcription-dependent neuronal plasticity: the nuclear calcium hypothesis. Eur. J. Biochem. 267, 5280–5283.

Pubmed Abstract | Pubmed Full Text

Balls, G. T., and Baden, S. B. (2004). “A large scale monte carlo simulator for cellular microphysiology,” in Proceedings of the 18th International Parallel and Distributed Processing Symposium (San Diego, CA: California University). doi: 10.1046/j.1432-1327.2000.01565.x

CrossRef Full Text

Bastian, P., Birken, K., Johannsen, K., Lang, S., Neuss, N., Rentz-Reichert, H., et al. (1997). UG a flexible software toolbox for solving partial differential equations. Comp. Vis. Sci. 1, 27–40. doi: 10.1007/s007910050003

CrossRef Full Text

Bentley, J. L. (1975). Multidimensional binary search trees used for associative searching. Commun. ACM 18, 509–517. doi: 10.1145/361002.361007

CrossRef Full Text

Blackwell, K. T., Wallace, L. J., Kim, B., Oliveira, R. F., and Koh, W. (2013). Modeling spatial aspects of intracellular dopamine signaling. Methods Mol. Biol. 964, 61–75. doi: 10.1007/978-1-62703-251-3_5

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Bower, J., and Beeman, D. (1997). The Book of GENESIS: Exploring Realistic Neural Models with the GEneral NEural SImulation System. New York, NY: Springer.

Brandi, M., Brocke, E., Talukdar, H. A., Hanke, M., Bhalla, U. S., Kotaleski, J. H., et al. (2011). Connecting MOOSE and NeuroRD through MUSIC: towards a communication framework for multi-scale modeling. BMC Neurosci. 12(Suppl. 1):P77. doi: 10.1186/1471-2202-12-s1-P77

CrossRef Full Text

Broser, P. J., Schulte, R., Roth, A., Helmchen, F., Waters, J., Lang, S., et al. (2004). Nonlinear aniostropic diffusion filtering of three-dimensional image data from 2-photon microscopy. J. Biomed. Opt. 9, 1253–1264. doi: 10.1117/1.1806832

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Burette, A. C., Lesperance, T., Crum, J., Martone, M., Volkmann, N., Ellisman, M. H., et al. (2012). Electron tomographic analysis of synaptic ultrastructure. J. Comp. Neurol. 520, 2697–2711. doi: 10.1002/cne.23067

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Cannon, R. C., O Donnel, C., and Nolan, M. F. (2010). Stochastic ion channel gating in dendritic neurons: morphology dependence and probabilistic synaptic activation of dendritic spikes. PLoS Comp. Biol. 6:e1000886. doi: 10.1371/journal.pcbi.1000886

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Carnevale, N. T., and Hines, M. L. (2006). The NEURON Book. Cambridge UK, Cambrdige University Press. doi: 10.1017/CBO9780511541612

CrossRef Full Text

Chen, X., Winters, C., Azzam, R., Li, X., Galbraith, J. A., Leapman, R. D., et al. (2008). Organization of the core structure of the postsynaptic density. Proc. Natl. Acad. Sci. U.S.A. 133, 4453–4458. doi: 10.1073/pnas.0800897105

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Clapham, D. E. (2007). Calcium signaling. Cell 131, 1047–1058. doi: 10.1016/j.cell.2007.11.028

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Cormen, T. H., Leiserson, C. E., and Rivest, R. L. (2011). Introduction to Algorithms. Cambridge, MA: MIT Press; Massachusetts Institute of Technology.

Deuflhard, P., and Weiser, M. (2011). Numerische Mathematik 3. Berlin: De Gruyter Lehrbuch 2011.

Djurfeldt, M., Hjorth, J., Eppler, J. M., Dudani, N., Helias, M., Potjans, T. C., et al. (2010). Run-time interoperability between neural network simulators based on the MUSIC framework. Neurinform 8, 43–60. doi: 10.1007/s12021-010-9064-z

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Egelman, D. M., and Montague, P. R. (1999). Calcium dynamics in the EC space of mammalian neural tissue. Biophys. J. 76, 1856–1867. doi: 10.1016/S0006-3495(99)77345-5

CrossRef Full Text

Eggermann, E., Bucurenciu, I., Goswami, S. P., and Jonas, P. (2011). Nanodomain coupling between Calcium channels and sensors of exocytosis at fast mammalian synapses. Nat. Rev. Neurosci. 13, 7–21. doi: 10.1038/nrn3125

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Eymard, R., Gallouet, T. R., and Herin, R. (2000). “Parabolic equations, Chapter 4, in The Finite Volume Method in Handbook Of Numerical Analysis,” Vol. 7. eds P. G. Ciarlet and J. L. Lions (Elsevier), 713–1020.

Gewaltig, M. O., and Diesmann, M. (2007). Nest (neural simulation tool). Scholarpedia 2, 1430.

Graham, L. (1999). Interpretations of Data and Mechanisms for Hippocampal Pyramidal Cell Models. New York, NY: Plenum Publishing Corporation.

Grillo, A., Logashenko, D., Stichel, S., and Wittum, G. (2010). Simulation of density-driven flow in fractured porous media. Adv. Water Resour. 33, 1494–1507. doi: 10.1016/j.advwatres.2010.08.004

CrossRef Full Text

Hackbusch, W. (1985). Multi-Grid Methods. Leipzig: Springer.

Hackbusch, W. (2003). Multi-Grid Methods and Applications in Springer Computational Mathematics. Leizpig: Springer.

Hansen, S., Henning, A., Ngel, A., Heisig, M., Wittum, G., Neumann, D., et al. (2008). In-silico model of skin penetration based on experimentally determined input parameters. Part i: Experimental determination of partition and diffusion coefficients. Eur. J. Pharm. Biopharm. 68, 352–367. doi: 10.1016/j.ejpb.2007.05.012

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Hayashi, Y., and Majewska, A. (2005). Dendritic spine geometry: functional implication and regulation. Neuron 46, 529–532. doi: 10.1016/j.neuron.2005.05.006

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Helmchen, F., Imoto, K., and Sakmann, B. (1996). Ca2+ buffering and action potential-evoked ca2+ signaling in dendrites of pyramidal neurons. Biophys. J. 70, 1069–1081. doi: 10.1016/S0006-3495(96)79653-4

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Heppner, I., Lampe, M., Nägel, A., Reiter, S., Rupp, M., Vogel, A., et al. (2013). “Software framework uG4: parallelmultigrid on the hermit supercomputer,” in High Performance Computing in Science and Engineering, eds W. E. Nagel, D. H. Kröner, and M. M. Resch (Berlin-Heidelberg: Springer), 434–449.

Hines, M. L., and Carnevale, N. T. (2003). “The NEURON simulation environment,” in The Handbook of Brain Theory and Neural Networks, ed M. A. Arbib Vol. 2, (Cambridge MA: MIT Press), 769–773.

Hodgkin, A. L., and Huxley, A. F. (1952). A quantitative description of membrane currents and its application to conduction and excitation in nerve. J. Physiol. 117, 550–544.

Pubmed Abstract | Pubmed Full Text

Ierusalimschy, R., Celes, W., and de Figueiredo, L. H. (2013). Lua: The Programming Language. Rio: MIT License.

Jungblut, D., Queisser, G., and Wittum, G. (2011). Inertia based filtering of high resolution images using a gpu cluster. Comp. Vis. Sci. 14, 181–186. doi: 10.1007/s00791-012-0171-2

CrossRef Full Text

Kanamori, T. (2013). Compartmentalized calcium transients trigger dendrite pruning in Drosophila sensory neurons. Science 340, 1475–1478. doi: 10.1126/science.1234879

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Katona, G., Kaszas, A., Turi, G. F., Hajos, N., Tamas, G., Vizi, E. S., et al. (2011). Roller coaster scanning reveals spontaneous triggering of dendritic spikes in CA1 interneurons. Proc. Natl. Acad. Sci. U.S.A. 108, 2148–2153. doi: 10.1073/pnas.1009270108

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Kerr, R. A., Bartol, T. M., Kaminsky, B., Dittrich, M., Chang, J. C., Baden, S. B., et al. (2008). Fast monte carlo simulation methods for biological reaction-diffusion system in solution and on surfaces. SIAM J. Sci. Comput. 30, 3127–3249. doi: 10.1137/070692017

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Kits, K. S., Dreijer, A. M. C., Lodder, J. C., Borgdorff, A., and Wadman, W. J. (1997). High intracellular calcium levels during and after electrical discharges in molluscan peptidergic neurons. Neuroscience 79, 275–284. doi: 10.1016/S0306-4522(96)00651-3

CrossRef Full Text

Korkotian, E., and Segal, M. (1999). Release of calcium from stores alters the morphology of dendritic spines in cultured hippocampal neurons. Proc. Natl. Acad. Sci. U.S.A. 96, 12213–12215. doi: 10.1073/pnas.96.21.12068

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Lee, D. T. (1997). Worst-case analysis for region and partial region searches in multidimensional binary search trees and balanced quad trees. Acta Inf. 9, 23–29.

Maravall, M., Mainen, Z. F., Sabatini, B. L., and Svoboda, K. (2000). Estimating intracellular calcium concentrations and buffering without wavelength ratioing. Biophys. J. 78, 2655–2667. doi: 10.1016/S0006-3495(00)76809-3

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

McDougal, R. A., Hines, M. L., and Lytton, W. W. (2013a). Reaction-diffusion in the NEURON simulator. Front. Neuroinform. 7:28. doi: 10.3389/fninf.2013.00028

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

McDougal, R. A., Hines, M. L., and Lytton, W. W. (2013b). Water-tight membranes from neuronal morphology files. J. Neurosci. Methods 220, 167–178. doi: 10.1016/j.jneumeth.2013.09.011

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Migliore, M., Morse, T. M., Davison, A. P., Marenco, L., Shepherd, G. M., and Hines, M. L. (2003). ModelDB: making models publicly accessible to support computational neuroscience. Neuroinformatics 1, 135–139. doi: 10.1385/NI:1:1:135

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Milner, B., L.R., S., and E.R., K. (1998). Cognitive neuroscience and the study of memory. Neuron 20, 445–468. doi: 10.1016/S0896-6273(00)80987-3

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Mount, D., and Arya, S. (2012). ANN Library.

Muha, I., Grillo, A., Heisig, M. Schöneberg, M., Linke, B., and Wittum, G. (2012). Mathematical modeling of process liquid flow and acetoclastic methanogenesis under mesophilic conditions in a two-phase biogas reactor. Biores. Technol. 106, 1–9. doi: 10.1016/j.biortech.2011.11.087

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Muha, I., Nägel, A., Stichel, S., Grillo, A., Heisig, M., and Wittum, G. (2011). Effective diffusivity in membranes with tetrakaidekahedral cells and implications for the permeability of human stratum corneum. J. Memb. Sci. 368, 18–25. doi: 10.1016/j.memsci.2010.10.020

CrossRef Full Text

Muller, D., Nikonenko, I., Jourdain, P., and Alberi, S. (2002). LTP, memory and structural plasticity. Curr. Mol. Med. 2, 605–611. doi: 10.2174/1566524023362041

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Murase, S., and Schuman, E. M. (1999). The role of cell adhesion molecules in synaptic plasticity and memory. Curr. Opin. Cell Biol. 11, 549–553. doi: 10.1016/S0955-0674(99)00019-8

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Nägel, A., Hansen, S., Neumann, D., Lehr, C. M., Schaefer, U. F., Wittum, G., et al. (2008). In-silico model of skin penetration based on experimentally determined input parameters. Part ii: Mathematical modelling of in-vitro diffusion experiments. identification of critical input parameters. Eur. J. Pharma. Biopharm. 68, 368–379. doi: 10.1016/j.ejpb.2008.11.009

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Nägel, A., Heisig, M., and Wittum, G. (2009). A comparison of two- and three-dimensional models for the simulation of the permeability of human stratum corneum. Eur. J. Pharm. Biopharm. 72, 332–338. doi: 10.1016/j.ejpb.2008.11.009

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Nagerl, U. V., Mody, I., and Vergara, J. L. (2000). Binding kinetics of calbindin-d(28k) determined by flash photolysis of caged ca(2+). Biophys. J. 79, 3009–3018. doi: 10.1016/S0006-3495(00)76537-4

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Oliveira, R., Terrin, A., Di Benedetto, G., Cannon, R. C., Koh, W., Kim, M., et al. (2010). The role of type 4 phosphodiesterases in generating microdomains of cAMP: large scale stochastic simulations. PLoS ONE 5:e11725. doi: 10.1371/journal.pone.0011725

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Popov, V. I., Kleschevnikov, A. M., Klimenko, O. A., Stewart, M. G., and Belichenko, P. V. (2011). Three-dimensional synaptic ultrastructure in the dentate gyrus and hippocampal area ca3 in the ts65dn mouse model of down syndrome. J. Comp. Neurol. 519, 1338–1354. doi: 10.1002/cne.22573

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Pumplin, D. W., Reese, T. S., and Llinas, R. (1981). Are the presynaptic membrane particles the calcium channels? Proc. Natl. Acad. Sci. U.S.A. 11, 7210–7213. doi: 10.1073/pnas.78.11.7210

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Queisser, G., Wittmann, M., Bading, H., and Wittum, G. (2008). Filtering, reconstruction, and measurement of the geometry of nuclei from hippocampal neurons based on confocal microscopy data. J. Biomed. Opt. 13:014009. doi: 10.1117/1.2829773

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Ray, S., and Bhalla, U. S. (2008). PyMOOSE: interoperable scripting in Python for MOOSE. Front. Neuroinf. 2:6. doi: 10.3389/neuro.11.006.2008

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Resasco, D. C., Gao, F., Morgan, F., Novak, I. L., Schaff, J. C., and Slepchenko, B. M. (2012). Virtual cell: computational tools for modeling in cell biology. Wiley Interdiscip. Rev. Syst. Biol. Med. 4, 129–140. doi: 10.1002/wsbm.165

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Shewchuk, J. R. (2002). “What is a good linear finite element? - interpolation, conditioning, anisotropy and quality measures,” Proceedings of the 11th International Meshing Roundtable (Berkeley, CA).

Si, H. (2009). TetGen: a quality tetrahedral mesh generator and three-dimensional delaunay triangulator. Int. J. Num. Methods Eng. 75, 856–880.

Silver, R. A., Lamb, A. G., and Bolsover, S. R. (1990). Calcium hotspots caused by l-channel clustering promote morphological changes in neuronal growth cones. Nature 343, 751–754. doi: 10.1038/343751a0

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Spacek, J., and Harris, K. M. (1997). Three-dimensional organization of smooth endoplasmic reticulum in hippocampal CA1 dendrites and dendritic spines of the immature and mature rat. J. Neurosci. 24, 4233–4241.

Pubmed Abstract | Pubmed Full Text

Tada, T., and Sheng, M. (2006). Molecular mechanisms of dendritic spine morphogenesis. Curr. Opin. Neurobiol. 16, 95–101. doi: 10.1016/j.conb.2005.12.001

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Tai, C., Kim, S., and Schuman, E. (2008). Cadherins and synaptic plasticity. Curr. Opin. Cell Biol. 5, 567–575. doi: 10.1016/j.ceb.2008.06.003

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Takahashi, A., Camacho, P., Lechleiter, D. J., and Herman, B. (1999). Measurement of intracellular calcium. Physiol. Rev. 79, 1089–1125.

Pubmed Abstract | Pubmed Full Text

Thompson, J. F., Soni, B. K., and Weatherill, N. P. (2012). Handbook of Grid Generation. Leipzig: CRC Press.

Van Aelst, L., and Cline, H. (2004). Rho GTPases and activity-dependent dendrite development. Curr. Opin. Neurobiol. 14, 297–304. doi: 10.1016/j.conb.2004.05.012

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Vogel, A., Reiter, S., Rupp, M., Nägel, A., and Wittum, G. (in press). uG 4 - A novel flexible software system for the simulation of PDE-based models on high performance computers. Comput. Vis. Sci.

Wald, I., and Hvran, V. (2006). “On building fast kd-trees for ray tracing and on doing that in o(NlogN),” IEEE Symposium on Interactive Ray Tracing (Amsterdam). doi: 10.1109/RT.2006.280216

CrossRef Full Text

Weinan, E., Engquist, B., and Huang, Z. (2003). Heterogenous multiscale method: a general methology for multiscale modeling. Phys. Rev. B 67, 367–450. doi: 10.1103/PhysRevB.67.092101

CrossRef Full Text

West, A., Griffith, E., and Greenberg, M. (2002). Regulation of transcription factors by neuronal activity. Nat. Rev. Neurosci. 3, 921–931. doi: 10.1038/nrn987

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Wikipedia (2012). K-d tree — Wikipedia, The Free Encyclopedia.

Wils, S., and De Shutter, E. (2009). STEPS: modeling and simulationg complex reaction-diffusion systems with python. Front Neuroinf. 3:15. doi: 10.3389/neuro.11.015.2009

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Wittmann, M., Queisser, G., Eder, A., Wiegert, J., Bengtson, C., Hellwig, A., et al. (2009). Synaptic activity induces dramatic changes in the geometry of the cell nucleus: interplay between nuclear structure, histone H3 phosphorylation, and nuclear calcium signaling. J. Neurosci. 29, 14687–14700. doi: 10.1523/JNEUROSCI.1160-09.2009

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Wolf, S., Grein, S., and Queisser, G. (2013). Employing NeuGen 2.0 to automatically generate realistic morphologies of hippocampal neurons and neural networks in 3D. Neuroinformatics 11, 137–148. doi: 10.1007/s12021-012-9170-1

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Xylouris, K. (2013). Generalized cable equation for signal processing in neurons. PhD Thesis.

Xylouris, K., Queisser, G., and Wittum, G. (2007). A three-dimensional mathematical model of active signal processing in axons. Comput. Visual. Sci. 13, 409–418. doi: 10.1007/s00791-011-0155-7

CrossRef Full Text