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

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

Ω 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

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:

as well as the conservation law for the buffer concentration:

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

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.

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).

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).

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

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

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.

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.

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

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.

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.

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.

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

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

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:

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:

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

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:

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

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

and

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

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

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:

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

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

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

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

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

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:

with boundary conditions

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

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