Multiscale modelling via split-step methods in neural firing

ABSTRACT Neuronal models based on the Hodgkin–Huxley equation form a fundamental framework in the field of computational neuroscience. While the neuronal state is often modelled deterministically, experimental recordings show stochastic fluctuations, presumably driven by molecular noise from the underlying microphysical conditions. In turn, the firing of individual neurons gives rise to an electric field in extracellular space, also thought to affect the firing pattern of nearby neurons. We develop a multiscale model which combines a stochastic ion channel gating process taking place on the neuronal membrane, together with the propagation of an action potential along the neuronal structure. We also devise a numerical method relying on a split-step strategy which effectively couples these two processes and we experimentally test the feasibility of this approach. We finally also explain how the approach can be extended with Maxwell’s equations to allow the potential to be propagated in extracellular space.


Introduction
Neurons are responsible for encoding information in the central nervous system.Lower-level functions are many times gathered at specific parts of the brain, processing information received from neurons throughout the body, often by means of signalling nerves, which are bundles of axons that reach out from the neurons.Higher-level functions depend on remarkably larger and complex networks of neurons with various types of feedback loops.
The chemical connections between neurons are handled via synapses, where neurotransmitters are extruded into the extracellular space by the presynaptic neuron.The neurotransmitters form a part of a chemical process which may initiate a potential wave into the postsynaptic neuronthis wave is known as the action potential.During the action potential, the membrane potential quickly rises and falls, and the resulting signal propagates along the cell [1].The process that underlies this propagation is the regulation of ion concentrations in both the intracellular cytoplasm and the extracellular space caused by integral membrane proteins called ion channels.
There are many variants of ion channel proteins, whose functions are only recently better understood through studies using experimental techniques such as X-ray crystallography [2].To our current understanding, ion channels reside at one of many conformal states, which are either closed (non-conducting) or open (conducting).If the conformal state is open, ions are allowed to pass through pathways called pores from the extrato intracellular space, or in the opposite direction.If the conformal state is closed, ions are blocked from entering the channel [3].Ion channels become activated either in response to a chemical ligand binding, or in response to voltage changes on the membrane [1,3], so-called voltage-gated ion channels.In this work, we focus on voltage-gated ion channels, which are important for the initiation and the propagation of action potentials along the neuronal fibre.
Mathematical models of neurons were initially formed by measuring electrically induced responses of a neuron using techniques such as the voltage or current clamp.From those experiments, transition properties and parameters of single channel gating could be identified [3].
Historically, models of firing neurons were formulated as ordinary differential equations (ODEs).However, in some specific cases experimental as well as theoretical findings suggest that the gating of ion channels is more accurately described by a stochastic process.The variance of the gating process, also known as the channel noise, is thought to be important for information processing in the dendrites and explain different phenomena regarding action potential initiation and propagation [4][5][6][7].For example, it has been shown that intrinsic noise is essential for the existence of subthreshold oscillations in stellate cells [8], that it contributes to irregular firing of cortical interneurons [9], and that it can explain firing correlation in auditory nerves [10].
Another important aspect of neuronal modelling is the investigation of action potentials propagating outward the neuron into the extracellular space.Simulation of such extracellular fields is one of the important methods used in computational neuroscience [11].Common usage includes the validation of experimental methods such as EEG and extracellular spike recordings, or in the modelling of physiological phenomena which can not be easily investigated empirically [12].
In this paper, we present a novel three-stage multiscale modelling framework consisting of the following components: (a) on the microscale, the gating process of the ion channels is governed by a continuous-time discrete-state Markov chain, (b) on the intermediate scale, the current-balance and the cable equation which are responsible for the action potential initiation and propagation is integrated in time as an ODE, (c) on the macroscale, the propagation of the trans-membrane current into an electrical field in extracellular space is achieved using partial differential equations (PDEs).
In x 2, the three modelling layers are explained in some detail.The numerical method via split-step modelling is summarized in x 3, where we also explain how the spatial coupling is achieved efficiently.We illustrate our method by some relevant examples in x 4 and offer a concluding discussion in x 5.

Neuronal modelling at different scales
We now describe the modelling framework at the individual scales, including the associated modelling assumptions.The microscale physics in the form of continuous-time Markov chains (CTMC) for the ion channels and the associated ion currents are discussed in x 2.1, the ODE-model for the action potential along the neuronal geometry in x 2.2, and the PDE-model of the extracellular electric field via Maxwell's equations in x 2.3.For convenience, a schematic explanation of the modelling framework is found in Figure 1 2

.1 Microscale: ion channel gating
The gating of voltage-dependent ion-channels is modelled by a Markov process.Hodgkin and Huxley [13,14] first proposed that the gating is governed by a set of gating variables.We assume that each gating variable takes values in a discrete state space and that a single combination of gating variables corresponds to an open and conducting ion channel [3].
The gating states and the transitions between them can be written in the form of a kinetic scheme.An example for the voltage-gated m 3 h 1 sodium channel is depicted in Figure 2 The notation for the channel states indicates that there are two involved gating variables m and h, where variable m takes one of four states and, independently, variable h takes one of two states.In total, we thus arrive at 8 states with a total of 10 reversible transitions, see Figure 2. The total number of channels in the open state m 3 h 1 implies a certain conductivity as will be detailed below.The transition rates between the states depend on the membrane voltage and on the state itself.When these transitions take place in a microscopic environment where molecular noise is present, a CTMC is the most suitable model.
We rewrite such schemes in a general, yet more manageable notation as follows.Assign the different combinations of gating variables an individual state S i , i ¼ 1; 2; . . .; M states .For the sodium scheme in Figure 2, there are 8 states and only the state S 8 corresponds to an open and conducting ion channel.The resulting equivalent scheme is shown in Figure 3 and an exemplary set of transitions and rates is (2:1) (2:2) where V m is the membrane potential at the location of the ion channel.For this particular example, the transition rates used by Hodgkin and Huxley are [13] 1 À e ÀðV m þ40Þ=10 ; (2:4) (2:5) (2:6) (2:7) The rates (2.4)-(2.7)were originally obtained by empirical fitting to experimental data and are formulated relative to the resting potential of the particular neuron model.We consider a small neural compartment a, obtained from a discretization of the neuronal fibre into small enough segments such that the potential is approximately constant within each compartment.Let there be a total of N a ion channels of the considered type Figure 3.An equivalent scheme to Figure 2, where each combination of gating variables corresponds to a discrete state of a Markov chain.(e.g. the m 3 h 1 sodium channel) on the compartmental surface.At any instant in time t, let s i be the number of channels in state S i , i ¼ 1; 2; . . .; M states .The scheme in Figure 3 now directly translates into Markovian transition rules for the states s 2 Z 8 þ , i.e. the ionchannel counts, thus defining a CTMC ðt; To model the electrophysiological properties of the ion channel, we define the single open channel conductance γ, as well as the neuronal membrane area A a such that the area density of ion channels is % a ¼ N a =A a .If O a is the count of the open ion channels, e.g. for the above example this is O a ¼ S a 8 , we arrive at the total conductance Next, the trans-membrane current produced by the ion channel is given by Ohm's law, where E is the reverse potential of the ion channel, that is, the membrane potential under which the trans-membrane current I a is zero, and where V a m is the membrane potential in compartment a.
In the macroscopic setting proposed by Hodgkin and Huxley [13], the transitions between the gating states happened according to first order reaction kinetics.The fraction of closed channels ð1 À o a Þ becomes open with rate α and the fraction of open channels o a becomes closed with rate β.Note that α and β may be voltagedependent for specific types of ion-channels, as in (2.1)-(2.7).The resulting ODE for each macroscopic gating variable o a can then be written as (2:10) Under the macroscopic formulation, the ratio of open channels is approximated as where M is the number of gating particles for the particular ion channel model, and p i is the exponent of the ith gating variable o i .As a concrete example, the sodium channel m 3 h 1 depicted in Figure 2 is gated by gating variables m and h, which in the macroscopic formulation have an exponent of three and one, respectively.In relation to (2.11), this means that o 1 ¼ m with p 1 ¼ 3, and o 2 ¼ h with p 2 ¼ 1, where m and h are here just symbols denoting the concentration of the gating variables.The ratio of open channels is then approximated as and the deterministic gating function for each variable obeys where α and β are the voltage dependent transition rates defined in (2.4)-(2.7).Solving (2.13)-(2.14)and using (2.9) thus constitutes the classical Hodgkin-Huxley ODE model for the ionic current.By rather evolving the Markov chain defined implicitly in Figure 3 and using (2.8)-(2.9),a stochastic current is defined.This stochastic model is the result of the microphysical assumption of discrete ion states obeying Markovian (memoryless) transition rules.For large-enough numbers of ion channels, a transition into the corresponding deterministic model can be expected under rather broad conditions [15].When a fine enough compartmentalization of a particular neuron is made, however, the stochastic model can be expected to be more realistic in the sense that the effects missing in a corresponding deterministic model can not be ignored.
We now proceed to discuss the appropriate model for evolving the resulting current, be it stochastic or deterministic, over the neuronal morphology.

Intermediate scale: current-balance and cable equation
On the intermediate scale, we solve for the current-balance equation of the neuron, describing the relation between the membrane voltage and the different current sources.We first take the ionic current sources from the microscale into account.We next model the capacitance of the membrane, a so-called passive neuronal property.As the neuronal morphology is divided into several compartments, we also include a term for the propagation of voltage between the compartments, resembling the cable equation.
As before, denote the compartments by the index a 2 f1; . . .; M compartments g.For each connection between compartments, an entry is made in the adjacency matrix A [16].Typically, a compartment is connected to two or three neighbours.If a compartment is connected to only a single neighbour, it represents one of the end points of the neuron.An example is found in Figure 4.The current flowing in each compartment is composed from different components.The current source of compartment a is the total axial current with respect to its connected compartments b 2 BðaÞ.Note that as A represents a directed graph, the connected component b 2 BðaÞ represents in-neighbours that are connected via an edge in the direction to vertex a, as well as out-neighbours which are connected by an edge outwards of vertex a.According to Ohm's law the current equals where G a b is the conductance between compartments a and b, and where V a m is the trans-membrane potential of compartment a.
We then add the ionic currents computed at the microscale.For the sake of generality, we extend the previous notation of a single type of ion channel into c types, T ¼ T 1 ; . . .; T c f g , acting simultaneously.Generalizing (2.9) we arrive at the total ionic current We now deviate slightly from the discussion and first consider a space-continuous solution obtained by solving for the propagation of the potential on a line oriented along the x-axis.Adding the axial current flux I axial , the total trans-membrane current flux I m , and a possibly additional external current source I inj gives where C m is the capacitance of the neuronal membrane.We have that the gradient of the inter-compartmental leakage current I L is The leakage current accounts for the ions that diffuse out of the intracellular space due to the natural permeability of the neuronal membrane.Given a suitable resting potential, Ohm's law is Given the parabolic character of this equation, the Crank-Nicolson scheme [17] is a natural choice when designing numerical methods.
Turning now to the compartment-discrete version, the leakage current in each compartment a is given as where the membrane leakage conductivity G a L and the leak reverse potential E a L are regarded as constants.
The total trans-membrane current flux now becomes In order to solve for the membrane potential in each compartment a, we start from a discrete version of (2.17), using (2.15) and (2.22), To summarize, on the intermediate scale we insert the stochastic conductivity G a T computed from (2.8) into (2.16) and solve the current-balance equation (2.23) for each of the neural compartments.We simultaneously solve for the trans-membrane current defined in (2.22), to be coupled into the macroscopic fieldnext to be discussed.

Macroscale: extracellular field potentials
To simulate spatially non-homogeneous distributions of electrical fields produced by single neurons and neuronal networks, we use an electrostatic formulation of Maxwell's equations to be discretized with finite elements.We solve for the electric field intensity E in terms of the electric scalar potential V, The relevant dynamic form of the continuity equation with current sources Q j is given by (2:26) with J and ρ the current density and electric charge density, respectively.Further constitutive relations include (2:27) and Ohm's law (2:28) in which D denotes the electric flux density.Finally, Gauss' law states that (2:29) Upon taking the divergence of (2.28) and using the continuity equation (2.26) we get (2:30) Rewriting the electric charge density using Gauss' law together with the constitutive relation (2.27) and finally applying the gauge condition (2.25) twice we arrive at the time-dependent potential formulation The values for the electric conductivity σ and the relative permittivity ε r were obtained from [18].The source currents Q j are computed from the compartmental model described in x 2.2.Specifically, in compartment a, we put where I a m ðtÞ is obtained by solving (2.22) and where we recall that A a is the area of the neuronal membrane in compartment a .
The boundary conditions here are homogeneous Neumann conditions (electric isolation) everywhere except in a single point which we take to be ground (V ¼ 0).In all our simulations this point was placed at the axis of rotation of the enclosing cylindrical extracellular space, and directly underneath the neuronal geometry.This procedure ensures that the formulation has a unique solution; without this specification it is otherwise only specified up to a constant.

Numerical modelling
The processes at the three scales, that is, the microphysics of the ion channel gating, the neuron current at the intermediate scale, and the propagation of the electric field, all take place in continuous time.Also, all processes formally affect each other through two-way couplings.
In x 3.1, we discuss the numerical coupling of the ion channel gating and the cable equation via a split-step strategy.In x 3.2, the propagation of the electric potential into extracellular space is summarized and here we also explain how the often rather complicated neuronal geometry is handled.To simplify matters, we will disregard the usually much weaker coupling from the external electric field back to the gating process.

Numerical coupling of firing processes
To get to grip with the details of the coupling between the ion channel gating process and the propagation of the action potential along the neuron, we need a more compact notation as follows.In compartment a, denote by X a ðtÞ the gating state, that is, the number of channels being in each of the different states at time t.For our sodium channel example, we have X a ðtÞ ¼ ½s a 1 ; s a 2 ; . . .; s a 8 ðtÞ T .The CTMC may be written compactly as dX a ðtÞ ¼ Sμ a ðX a ; V a m ; dtÞ; (3:1) for a = 1,. .., M compartments , and where S is an M states × N transitions matrix of integer transition coefficients.The dependency on the state (X a ;V a m ) is implicit in the random counting measure µ a associated with an N transitions -dimensional Poisson process.As a concrete example, the transition S3 → S4 in (2.1) satisfies with that is, with the understanding that (2.1) is the i th transition according to some given ordering.Since the process is a Poisson process, (3.2) expresses independent exponentially distributed waiting times of intensity α m ðV a m Þ for each of the X a 3 possible transitions of type S 3 !S 4 .
The propagation of the action potential can similarly be written in more compact ODE notation, where G is just the current-balance equation (2.23) and which depends on the state ðX a ; V a m Þ via (2.8), (2.16), and (2.22), respectively.Eqs.(3.1) and (3.4) form a coupled CTMC-ODE model which falls under the scope of Piecewise Deterministic Markov Processes (PDMPs) and for which numerical methods have been investigated [19].A highly accurate implementation is possible through event-detection in traditional ODE-solvers.This, however, has a performance drawback since the ODE-solver must continuously determine to sufficient accuracy what microevent happens when.In the experiments in x 4 we have chosen to rely on a simple splitstep strategy as follows.Given a discrete time-step Δt, and Sμ a ðX a ðsÞ; V a m;n ; dsÞ; GðX a nþ1 ; V a m ðsÞ; fV b m ðsÞg b2BðaÞ Þds: That is, (3.5) evolves the CTMC (3.1) keeping the voltage potential fixed at its value from the previous time-step t n .Similarly, (3.6) evolves the ODE (3.4) while keeping the state of the ion channels fixed at the time-step t nþ1 .Importantly, the usually more expensive stochastic simulation in (3.5) is fully decoupled as it only depends on the state of the compartment a.The global step is achieved separately by solving the connected ODE in (3.6), and is usually quite fast.The approximation (3.5)-(3.6)can be understood as a split-step method and may also be analysed as such [20,21].The order of the approximation can then be expected to be 1=2 in the root mean-square sense, Although it appears to be difficult to increase the stochastic order of this approximation, the accuracy is likely going to increase by turning to symmetric Strang-type splitting methods for the stochastic part [21], possibly also adopting a higher order scheme for the ODE-part.However, the efficiency of such an approach ultimately depends on the strength of the non-linear feedback terms and is difficult to analyse a priori.In the present proof-of-concept context, we aim at a convergent and consistent coupling for which (3.5)-(3.6)are a suitable starting point.We thus postpone the investigation of more advanced integration methods for another occasion.

Spatial extension of the firing process
The modelling in x 2.2 did not take into account a possible free-space potential V ext , external to the cell.The necessary modifications to incorporate this are as follows.For x 2 Ω & R 3 , let V ext ðxÞ be a given external potential and denote by V a ext the value at compartment a .Then replace (2.23) with for G a m the membrane transneuronal conductivity.With this modification, (2.23) propagates the effect of the given external field V ext ðxÞ along the neuron.
In x 2.3, we discussed how the effect of a trans-membrane current propagates into extracellular space as an electric potential V (cf.(2.31)-(2.32)).With the previous discussion in mind, the following two-way coupling thus emerges: the solution of the current-balance equation (2.23) yields a current source, which feeds into (2.32).In turn, this implies an external potential V ext : ¼ V, obtained by solving (2.31), finally to be inserted into (3.8)above.The full model coupling thus arrived at takes the schematic form It is certainly possible to devise a numerical method using a similar split-step strategy as in x 3.1 for this full coupling.As mentioned before, to simplify we shall assume that the feedback from the external field and onto the neuron is weak, such that V ext % 0 in (3.8) and we thus consider the simplified model (see also Figure 1) This assumption is valid whenever the simulation consists of a small number of neighbouring neurons, but becomes inaccurate if a large number of nearly parallel neurons is considered.See [22,23] for a further discussion.It follows that the solution of the field potential V can be done 'offline'.That is, this problem may be solved in isolation using a pre-recorded current source I a m ðtÞ obtained from the split-step method described in x 3.1.
The three-dimensional neuronal geometry was constructed in Comsol Multiphysics with the help of the interface to Matlab ('LiveLink') by morphological additions and Boolean unification of simple geometric objects, representing neuronal compartments.The TREES Toolbox [16] was used to conveniently access the geometrical properties of single compartments over the Matlab interface.More specifically, the geometry was constructed by parsing the adjacency graph A .Starting from the initial node of the graph, this procedure makes sure that each compartment is connected to the same neighbours as in the numerical model for the axial current flow in (2.15).
In an initial attempt, we aimed to represent the 3D geometry as an exact counterpart of the compartmental model, where each compartment is understood as a cylinder with a certain length and diameter.If the direction of the main axes of any two joining cylinders differed, a sphere was added in between the cylinders, followed by a removal of all interior boundaries.Although this created a direct volumetric representation of the neuronal compartments, the approach is difficult to generalize to neuronal branches with a more complicated connectivity.The reason is that the triangulation of the final object becomes extremely difficult to achieve as the mesh engine insists on fully resolving the curvature.See Figure 5 for an example of a problematic mesh, emerging at the intersection of a cylinder and a sphere.
A second and more successful attempt was made by which three-dimensional curves made up of line segments for each compartment was constructed.This simplifies the meshing process immensely, since the extracellular mesh is not constrained by highcurvature cylindrical boundaries.Implicit here is the assumption that the neuron is very thin compared to the external length scale of practical interest.This is valid as the diameter of a dendritic structure is % 1m and the considered extracellular space is typically in the range of several mm [11].In Figure 6 we show examples of both approaches.

Experiments
We devote this section to some feasibility experiments of the method proposed.In x 4.1, we look in some detail at the coupling of the microscopic and intermediate scales.
Using the stochastic currents so computed, the induced electrical field is propagated outside the neuron in x 4.2.

Micro-meso coupling
In this section, we provide numerical experiments of the coupling between the microscopic scale, introduced in x 2.1, and the mesoscopic scale, introduced in x 2.2.We simulate the classical squid model proposed by Hodgkin and Huxley [13], which is also a part of a widely used benchmark for neuronal simulators [24].The model includes two types of voltage-gated ion channels, Na þ and K þ ions, which in our case are both modelled as continuous-time Markov processes (cf.x 2.1).The kinetic gating scheme for the Na þ channel is shown in Figure 2. The K þ channel follows a similar scheme, but has only four discrete gating states S 1 ; . . .; S 4 , of which only the state S 4 is open.The transition rates in this case are voltage dependent as follows [13] Note the similar form of (4.1)-(2.4)and (4.2) -(2.5), respectively, which is due to using the same empirical fitting procedure.The density of the ion channels is % K ¼ 30 m À2 and % Na ¼ 330 m À2 , respectively [3].The single channel conductance equals γ K ¼ 360=30 pS, and γ Na ¼ 1200=330 pS, while the reversal potentials are E K ¼ À77 mV and E Na ¼ 50 mV [24].
The parameters for the intermediate scale are as follows.The specific membrane capacitance c m ¼ 1FðcmÞ À2 , resting potential E r ¼ À65 mV, cytoplasm resistivity % c ¼ 100 Ωcm, specific leak conductance g L ¼ 40000 À1 SðcmÞ À2 and leak reversal potential E L ¼ À65 mV [24].We let the model be confined to a cylindrical geometry with a length of 1 mm and a diameter of 1 m.The cylinder is compartmentalized into M compartments sub-cylinders of equal length and diameter.The root node is ignited by a current injection of 0:1 nA.
In practice, we use the Gillespie's 'Direct Method' [25] to evolve (3.5) and the Crank-Nicholson scheme [17] to solve (3.6), relying on the fact that the latter is linear in V a m .In Figure 7, we show the numerical solution of the coupled model, overlaid with the deterministic solution where ion channels are described by ODEs as in (2.13)-(2.14).We show three traces of the membrane voltage V m in one neuronal compartment over time, where the injected current has been varied from a lower value (0:045nA) to a larger value (0:1nA).The dynamics of the stochastic model for the initial current injections clearly differ from the dynamics of the deterministic one, where a single spike or a train of spikes can be triggered in the stochastic representation, while no spike can be obtained in the deterministic one.At the higher amount of current injection, we observe a trace with similar characteristics for both model representations, but with an increasingly different phase shift.
In Figure 8, we show a numerical convergence study of the coupled model, concerning two method parameters.We show how the interspike interval (ISI) changes as a function of the coupling time step Δt, as well as the discretization of the geometry Δx.The ISI is defined as the duration between the peaks of two spikes.As the neuronal firing is now a stochastic process, the ISI will be a distribution, and hence we present the first and second moments.We find that, for the study of the coupling time step, the ISI appears to be well resolved at Δt ¼ 10 À1 ms for the spatial discretization presented.For the spatial discretization, we find that the ISI distributions do not significantly differ for a voxel length of under 1m.

Three-scale coupling
For this example, we took the morphological description of a pyramidal CA1 neuron from [26], which contains about 1500 compartments.We re-sampled the geometry in order to aggregate short compartments of sizes less than 50 m into larger compartments, leading to a final representation consisting of approximately 400 compartments.
Although the reduction might alter the properties of the model and more evolved protocols for compartment reduction could be used [27,28], we ignore the induced discretization error here since the purpose of the model is to demonstrate the overall numerical method.We solved (3.5) and (3.6) over the morphology, incorporating the active properties described in x 4.1.Next, we scaled the transmembrane currents I m ðtÞ (cf.(2.32)) and mapped this to the corresponding curve segment as a current source Q j ðtÞ [29].It can be noted in passing that we here assume that transmembrane currents are the only cause of change of extracellular potential, which is not the case in a real neuron, as for example synaptic calcium-mediated currents are suspected to contribute to a large fraction of the extracellular signature.
Equipped with the source currents (2.32), the formulation (2.31) is efficiently solved by Comsol's 'Time discrete solver', which is based on the observation that the variable W: ¼ ΔV satisfies a simple ODE.Solving for W in an independent manner up to some time t, it is then straightforward to solve a single static PDE to arrive at the potential V itself.For the Time discrete solver, the time-step was set to Δt as used in the split-step method, thus ensuring a correct transition to the macroscopic scale.A tetrahedral mesh [23] was applied to discretize space (using the 'finest' mesh setting; resolution of curvature 0.2, resolution of narrow regions 0.8).The simulations were verified against coarser mesh settings in order to ensure a practically converged solution.
The result of the simulation is visualized in Figure 9.We have inserted four 'point probes' at a radial distance of 1000 m to the neuron, measuring the extracellular voltage at single points.The electric potential thus monitored by the probes is shown in Figure 10.

Conclusions
In this paper, we have proposed a model framework for neuronal firing processes consisting of three layers: on the microscale, ion channel gating is modelled as a CTMC.On the intermediate scale, the currents produced by open channels are integrated into the current-balance equation proposed by Hodgkin and Huxley.Finally, on the macroscale, the outward current of neurons is propagated into extracellular space, simulating the emission of an extracellular potential.We have also described a numerical approach to the coupling of the different scales and indicated through computational results the feasibility of the overall approach.
To date, several exact and approximate simulation methods [4,[30][31][32] of stochastically gated ion channel models have been proposed.However, to our knowledge, the coupling between the microscopic gating layer and the mesoscopic layer describing the action potential initiation and propagation, has not yet been rigorously studied.In this paper, we provide the formulation of the coupling as a split-step method and moreover numerically observe the convergence of the method with respect to the size of the coupling time-step as well as the spatial discretization.We show that the distribution of interspike intervals may be significantly affected by the choice of the coupling time-step, but does not appear to depend strongly on the spatial discretization.Finally, we discuss the theory and praxis of incorporating the macroscopic scale into the model, including the spatial representation of neuronal compartments and appropriate numerical procedures for the simulation of the local field potential (LFP) propagation.
While it has often been taken for granted that ion channel gating occurs deterministically, accumulating research evidence indicates that the presence of stochasticity significantly influences the neuronal behaviour [5][6][7].In turn, neurons respond with a high level of variability to the repeated presentation of equal stimuli.This leads to remarkable difficulties in studying the link between single cell biophysical properties and their function in larger neuronal networks, both in health and disease.Additionally, a recent study [4] has shown that stochastic ion channel gating largely differs not only between different neuronal cell types, but also locally between different parts of the neuron.Given the technical difficulties of assessing the signal channel properties of smaller neuronal compartments such as axons and dendrites, developing reliable mathematical models is essential to tackle this problem.
Several studies [4][5][6][7]9,10] have demonstrated how incorporation of channel noise into the Hodgkin and Huxley equations could resemble multiple realistic neuronal behaviours.Although it is not the scope of this work to mimic any particular experimental problem, but rather to provide a modelling framework that could be used by a diverse group of scientists in the future, we foresee multiple interesting applications.For example, it was shown previously [7,33,34] that the stochastic nature of voltagegated ion channels in the medial entorhinal cortex stellate cells are crucial for generating subthreshold oscillations in theta ( % 8Hz) frequency range.These intrinsic oscillations in stellate cells were suggested to generate theta oscillations in vivo LFP [34].It is well established that theta oscillatory activity in vivo provides a temporal window in which spatial and declarative memories are formed [35].Thus, our framework could be used to test the link between the stochastic nature of single channels in specific cell types and their effect at both the cellular and the network level.

Figure 1 .
Figure 1.Schematic overview of the proposed multiscale modelling framework.The CTMC solution on the microscale depends on the membrane voltage that is computed on the mesoscale.This coupling is bidirectional; the mesoscale solution depends on the microscopic stochastic ion current.The macroscopic solution considered here depends solely on the mesoscopic membrane current.

Figure 2 .
Figure 2. Kinetic scheme for the m 3 h 1 sodium channel gating proposed by Hodgkin and Huxley [13].Only the gating state m 3 h 1 represents an open ion channel, for all other states the ion channel is closed.

Figure 4 .
Figure 4. Geometry of a neuron and the associated adjacency matrix A (figure adapted from [16]).A nonzero entry in A implies a connection between the corresponding compartments.

Figure 5 .
Figure 5. Example of a problematic mesh at the intersection of two cylindrical compartments.

Figure 6 .
Figure 6.Example of a neuronal morphology created with cylindrical objects (left), and with curves (right).

Figure 7 .
Figure 7. Membrane voltage behaviour around threshold current.The accuracy of the deterministic solver was verified with the reference solution of Rallpack 3. In all cases the discretization is Δt ¼ 0:05 and Δx ¼ 0:05 .Upper: injected current is I inj ¼ 0:045 nA, middle: injected current is I inj ¼ 0:063 nA, lower: injected current is I inj ¼ 0:1 nA.

Figure 8 .
Figure 8. Top: The mean interspike interval (ISI) AE standard error as a function of the time step Δt .The compartment length is here Δx ¼ 0:02 .Bottom: The ISI as a function of the compartment length Δx at a time step Δt ¼ 0:05 .The red line represents the interspike interval from the Rallpack 3 reference solution[24].The number of trajectories in all runs is N ¼ 40 .

Figure 9 .Figure 10 .
Figure 9. Solution of the full three-scale model framework: propagation of an extracellular action potential into homogeneous extracellular space.