Brain-wave equation incorporating axodendritic connectivity

We introduce an integral model of a two-dimensional neural ﬁeld that includes a third dimension representing space along a dendritic tree that can incorporate realistic patterns of axodendritic connectivity. For natural choices of this connectivity we show how to construct an equivalent brain-wave partial differential equation that allows for efﬁcient numerical simulation of the model. This is used to highlight the effects that passive dendritic properties can have on the speed and shape of large scale traveling cortical waves.


I. INTRODUCTION
Ever since Hans Berger made the first recording of the human electroencephalogram (EEG) in 1924 there has been a tremendous interest in understanding the physiological basis of brain rhythms [1]. This has included the development of so-called neural field models, recently reviewed in [2], for the forward generation of EEG signals as well as the use of techniques from spatiotemporal pattern formation [3], and nonequilibrium phase transitions [4] for their analysis. These neural field models are often formulated in terms of delayed integrodifferential equations posed on appropriate surfaces, including lines [5], planes [6], and folded cortical structures [7].
At heart modern biophysical theories assert that EEG signals from a single scalp electrode arise from the coordinated activity of ∼10 8 pyramidal cells in cortex [8]. These are arranged with their dendrites roughly in parallel and perpendicular to the cortical surface. Synaptic activation at the dendrites creates a net ionic membrane current representing a sink (source) for excitatory (inhibitory) synapses with a negative (positive) extracellular potential. Because there is no longterm accumulation of charge in the tissue, this synaptic current flowing inside the cell escapes across the membrane again as a return current, in particular in places with large surface area like the soma. This creates a distributed extracellular source in the case of a synaptic sink and vice versa for a synapse that acts as a source. Due to the elongated morphology of * james.ross1@nottingham.ac.uk † m_margetts@btinternet.com ‡ i.bojak@reading.ac.uk § rachel.nicks@nottingham.ac.uk ¶ d.avitabile@vu.nl ** stephen.coombes@nottingham.ac.uk pyramidal neurons, the separation of sink and source typically leads to the effective formation of a current dipole. Hence, at the population level the potential field generated by a synchronously activated population of cortical pyramidal cells behaves like that of a current dipole layer. Although the important contribution that single dendritic trees make to generating extracellular electric field potentials has been realized for some time, and can be calculated using Maxwell's equations [9], they are typically not accounted for in neural field models. For example in the Nunez model [10] for the generation of the α-rhythm neuronal cell types are arranged in sheets, with no representation of dendritic processing. A recent study of standing and traveling waves in the Nunez model posed on a sphere can be found in [11]. Interestingly the Nunez model has a representation of space-dependent delays arising from the finite speed of action potential propagation along axons, and for certain patterns of decaying strength of nonlocal connectivity it can be formulated as a damped inhomogeneous wave equation [12]. This brain-wave equation, and variants thereof, has played a major role in the interpretation of EEG signals since the 1970s [13][14][15][16][17]. Moreover, the local nature of such models means that they are amenable to analysis with standard numerical techniques for partial differential equations (PDEs), circumventing the challenges of evolving delayed integrodifferential models [12,18].
With the advent of laminar electrodes to record from different layers of cortex [19] it is now timely to build on the original work of Bressloff, reviewed in [20], and develop neural field models that incorporate a notion of dendritic depth. Moreover, given the benefits of the brain-wave equation it is sensible to look for generalizations that can incorporate both axons and dendrites, as well as the patterns of connectivity between them. This is the topic that we address in this paper. We show how to develop the cable modeling approach of Rall [21] to describe a firing rate cortical tissue model with axodendritic patterns of synaptic connectivity, allowing for a mixture of excitatory and inhibitory synapses across the model of the dendritic tree. In doing so we obtain a natural extension of two-dimensional neural field models to include a third dimension representing position along a dendritic cable. The firing rate in the somatic (cell body) layer is taken to be a smooth sigmoidal function of the cable voltage at the soma, which is in turn determined by the spatiotemporal pattern of synaptic currents on the cable.
In Sec. II we describe the formulation of the model in terms of a generalized neural field, prescribed by a given pattern of axodendritic connectivity. This tissue model is a continuum description of cell bodies arranged in a surface with each point fibrated by a simple one-dimensional model of a dendritic tree. The dynamics of the cell bodies is described with a delayed integrodifferential equation. In turn this is used to drive activity in the dendrite along connected structural axonal pathways. The dynamics of the dendritic tree is considered to be passive and is modeled using the PDE approach of cable theory. Next, in Sec. III, we exploit the fact that the strength of connections in large scale cortical structures is known to fall off exponentially with distance to reformulate the model as a set of coupled PDEs. These generalize the traditional brain-wave equation to include a notion of dendritic depth. Numerical simulations of this local PDE formulation are presented in Sec. IV, with our numerical scheme validated against an exact traveling front solution that can be constructed in the limit of a steep sigmoidal firing rate. These simulations are used to highlight the effect that known correlations between the distance of connection along the dendrite and separation between cell bodies can have on the speed and shape of traveling waves. Finally in Sec. V we discuss natural extensions of the work in this paper.

II. THE MODEL
We take as our model dendrite the standard onedimensional unbranched uniform cable equation, and distinguish between pyramidal cells and interneurons by introducing the labels P (for pyramidal) and I (for interneuron). The voltage V a (x, r, t ) at position x ∈ R along an infinite passive cable of neuron type a ∈ {P, I} with somatic coordinate r ∈ R 2 can then be written Here, τ a is recognized as the membrane time-constant of the dendrite, and D a as the cable diffusion coefficient. The electrotonic length √ D a τ a of a dendrite is typically in the range 0.1-1 mm. Here, I a (x, r, t ) is the synaptic input, which we shall split into an excitatory and inhibitory part: where V + = V + (x) and V − = V − (x) are positive and negative synaptic reversal potentials, respectively, that can vary across the dendritic cable. The conductance changes g ab evolve according to a slightly modified neural field prescription as [2,22] Here, η ab (t ) is an α-function synaptic filter where H is the Heaviside step function. The time-to-peak for the α-function η ab (t ) is α −1 ab , and the parameter α ab can be used to control the speed of the synapse. A fast synapse would have a typical time-to-peak of around 1ms, whilst a slow synapse would have a time-scale of around 100 ms. The function W ab (x, r, r ) describes the axodendritic connectivity pattern between populations a and b, while v ab is the velocity of the action potential propagated from population b to a. This speed can range from around 0.5 m/s in unmyelinated axons to 150 m/s in myelinated axons (in the peripheral nervous system), and typical values for cortico-cortical axonal speeds in humans are distributed, and appear to peak in the 5-10 m/s range [23] though speeds in callosal fibres can range from 7-19 m/s [18]. The field h a is taken as a measure of the somatic activity in population a, while f a describes the firing rate of population a. This latter function is often chosen to have a simple sigmoidal form, and we shall work with controls the steepness of f a , and θ a is a threshold. As a model of h a we shall take the somatic potential, namely the voltage on the cable at the point where the cell body lies. We fix this to be the coordinate where x = 0 so that h a (r, t ) = V a (0, r, t ). A schematic of the neural field model incorporating dendrites is shown in Fig. 1.
Assuming vanishing initial data, we may write the solution to Eq. (1) in the form where G a (x, t ) is the Green's function and the operator ⊗ denotes spatiotemporal convolution over the (x, t ) coordinates. Unfortunately Eq. (5) does not provide an explicit formula for V a (and hence h a ) since I a itself depends on V a . However, by repeated substitution and truncating terms at second order in the conductances we have that which shows that the dendrite acts to mix conductance changes (if input currents are shunted). The formula for h a in Eq. (7) is in fact the first two terms in a Neumann series expansion and see [20] for a further discussion of this in the context of dendritic modeling.
To model cortical anatomy in a biologically plausible fash- where r = |r|, d ab , κ ab ∈ R and W 0 ab > 0, which incorporates the fact that synapses are located further away from the soma as the separation between neurons increases [24]. For κ ab = 0 there is no correlation between somatic and dendritic coordinates and synapses occur at a fixed distance d ab from the soma. To complete the description of the neural tissue model we need only specify the form of the anatomical connectivity function w ab (r) that describes how the strength of the interactions changes with separation between cell bodies. Most long-range synaptic interactions are excitatory, with excitatory pyramidal cells sending their myelinated axons to other parts of the cortex. Inhibitory interactions, on the other hand, tend to be much more short-ranged. For excitatory connections between cortical areas in macaque monkeys the weight of connection between two areas decays approximately exponentially with their wiring distance, with a characteristic distance of ∼11 mm [25], and for a nice overview of brain structure and dynamics across scales we recommend the article by Wang and Kennedy [26]. Hence, it is natural to choose an exponential form: with σ aP > σ aI to respect the fact that in neocortex the extent of excitatory connections w aP is broader than that of inhibitory connections w aI .

III. GENERALIZED BRAIN-WAVE EQUATION
To numerically analyze the nonlinear integrodifferential equation model with space-dependent delays given in Sec. II, it is convenient to consider the development of a PDE model where evolution is expressed in terms of differential, rather than integral, operators. This is a common approach in analyzing neural field models that lack dendrites [27] because such models are numerically easier to simulate than non-local integral equations. This reformulation in terms of a PDE essentially relies on the specific choice of the kernel (9) having certain properties, and the main one being that its Fourier transform can be well approximated by a rational function.
We make use of the fact that an α-function is the Green's function of a linear differential operator to write Hence from Eq. (3) we obtain a PDE for g ab as Q ab g ab = ψ ab (10) with g ab (x, r, 0) = 0, ∂g ab /∂t| t=0 = 0 and where . Introducing a Fourier transform over the dendritic spatial coordinate allows us to represent a function φ(x, r, t ) in the form Using Eq. (8) we find that with Taking further Fourier transforms with respect to (r, t ) and introducing spectral parameters (k, ω) we exploit the convolution structure of Eq. (13) to obtain where ρ b = f b • h b and the Fourier transform of H ab may be evaluated as and A ab (p, ω) = 1/σ ab + iω/v ab + ipκ ab . After cross multiplying and using the long-wavelength approximation [namely, expanding H ab (p, k, ω) around k = 0] and taking the inverse Fourier transform we obtain a PDE in three spatial dimensions (two for the somatic sheet and one for the dendritic cable) where We interpret Eq. (17), coupled to Eqs. (10) and (7), as the natural generalization of previous brain-wave equations of Nunez type, such as those in [28][29][30], to include dendritic processing. To work without the small conductance assumption one should use the PDE (1), instead of Eq. (7) with h a (r, t ) = V a (0, r, t ). A similar analysis can be carried out for a model with a one-dimensional somatic space. Using z as the coordinate in the somatic space and choosing an anatomical connectivity function w ab (z) = exp(−|z|/σ ab )/(2σ ab ) we obtain, following the same approach as above, the PDE in two spatial dimensions (one for the somatic space and one for the dendritic cable): Unlike its higher dimensional counterpart this equation is exact since the long-wavelength approximation is not needed.

IV. NUMERICAL SIMULATIONS
To numerically evolve the new type of brain-wave equation described in Sec. III it is natural to consider a finitedifference scheme to approximate spatial derivatives using central difference operations for second order derivatives, backward difference operations for first order derivatives and time-stepping that can be performed using an adaptive solver. For the latter we use a routine from the JULIA package-DifferentialEquations.jl [31,32]. The spatial domain is discretized into uniform rectangles with size x × z for the two dimensional model and cuboids with size x × r 1 × r 2 for the three-dimensional model. The Fourier transforms in the derivation of the brain-wave equation rely on infinite spatial domains and therefore we need to consider carefully the boundary conditions employed in a practical setting. It is natural to consider a periodic somatic domain with a large period, using a ring for the model with a one-dimensional somatic space and a torus when the somatic space is twodimensional. The cable equation may be imposed on a large finite domain with x ∈ [X − , X + ], and it is natural to choose closed end boundary conditions (so that the current is zero), that can be implemented as Neumann boundaries at x = X ± . The Dirac-δ function in Eq. (19) can be approximated with a thin narrow Gaussian δ ε (x) = exp(−x 2 /ε 2 )/ √ πε 2 with a fixed value of ε as the small dendritic spatial discretization size x. Further discussion of the numerical method is provided in Appendix A.
To validate our numerical approach it is useful to consider a restriction of the model that allows for an exact solution, so that theory and numerics can be easily compared. One such reduction is to focus on the model in two dimensions given by Eq. (19), and treat only a single excitatory population with a Heaviside firing rate and with a large spatially uniform synaptic reversal potential. This latter restriction means that inputs are no longer state dependent, and, after absorbing a factor of V + within g PP , that we may write I P = g PP . In this case, and fixing κ PP = 0, it is possible to construct a traveling front solution whose speed c PP is given by the implicit solution of Here, θ P is the threshold in firing rate function f P (h) = H (h − θ P ), σ PP is the spatial length scale for the exponential decay of connection strength, v PP is the speed of the action potential, and d PP is the contact distance of the synapse on the dendrite (as measured from the soma). The functions G(·, λ) and η(λ) are the Laplace transforms t → λ of G P and η PP given, respectively, by Eqs. (6) and (4). These can be calculated explicitly as The derivation of Eq. (20) is given in Appendix B. For the following numerical simulations we remove physical units by measuring space in units of σ PP (the spatial scale for anatomical connectivity) and time in units of τ P (the membrane time-constant of the cable). Thus we replace D P = D · σ 2 PP /τ P , α PP = α/τ P , d PP = d · σ PP , and v PP = v · σ PP /τ P in our equations and re-scale x → x · σ PP , r → r · σ PP , and t → t · τ P . This means also derived variables like x are to be considered as rescaled in the following, e.g., x = 0.1 means 0.1 · σ PP . For concreteness we shall consider σ PP = 10 mm and τ P = 10 ms. Conveniently, in this case σ PP /τ P = 1 m/s and thus front wave speed c PP = c · σ PP /τ P has the same numeric value in units of m/s as the unit-free c. In figure captions we will show the physical units in brackets "(m/s)" as equivalent for our specific choice. Unless otherwise stated, we shall also fix the electrotonic distance √ D P τ P = √ Dσ PP to be 0.1σ PP = 1 mm, i.e., D = 0.01.
In Fig. 2 we show a plot of the theoretical wave front speed, using Eq. (20), against direct numerical simulations obtained using the numerical scheme described above. The periodic somatic spatial domain is chosen sufficiently large to allow for the full formation of fronts, whilst the size of the dendritic domain is chosen in relation to its electrotonic length such that |d − X ± | > √ D. There is excellent agreement between theory and simulations under parameter variation in both d (the dendritic contact distance of the synapse from the soma) and θ (the threshold in the Heaviside firing rate function). We take this as validation of the proposed finite difference scheme, and reuse the same numerical methodology to explore the behavior of the model in more general settings, with a sigmoidal firing rate, finite synaptic reversal potential, and non-zero correlation parameter κ.
In Fig. 3 we further consider the two-dimensional model (with a one-dimensional somatic space) using a sigmoidal firing rate. Here, we observe traveling waves in both the conductance and voltage. Activity clearly spreads along the somatic space and decays in the dendritic space (as expected from the diffusive nature of the cable equation). The effect of increasing κ on the speed of the front for a Heaviside firing rate is quantified in Fig. 4. The results from these numerical experiments show that the wave speed decreases with increasing κ (for all values of d). In Fig. 5 the effects of changing the time-to-peak of the synapse α −1 is shown. Here we see that the wave front speed decreases as the synapse becomes slower. In Fig. 6 we show that an increase in the excitatory synaptic reversal potential causes an increase in wave front speed. When considering the variation of the correlation strength κ we find that for κ = 0 the peak of the wave response along the dendrite occurs at a distance d from the soma, as expected. However, with an increase in κ this peak response shifts slightly away from d to a larger distance from the soma. We illustrate this effect in Fig. 7, where we also see an increase in the width of the activity that spreads into the dendrite away from the contact point at x = d. Finally, in Fig. 8 we consider the three-dimensional model (with a two-dimensional somatic space) using a sigmoidal firing rate. Here we show a traveling wave propagating from the center of the domain. Similarly as for the two-dimensional model, see Fig. 3, the activity clearly spreads through the somatic layer, and decays along the dendrite at large distances from the cell body.

V. DISCUSSION
In this paper we have introduced a generalization of the Nunez brain-wave equation to account for passive dendritic processing. This is potentially important for the understanding of real EEG data given that this is known to arise from the addition of dipole moments generated by sinks and sources on dendrites. The generalized model bridges multiple scales, combining local models of dendrites with global models of synaptic and firing rate activity with both distributed (synaptic) and space-dependent (axonal) delays, emphasizing that dendritic response (Green's) functions can lead to another form of distributed delay that shapes emerging spatiotemporal network patterns. As well as being relevant to EEG the finer detail of the local model (with a spatially extended dendrite) makes it relevant to recorded potentials arising in electrocorticography (ECoG) and local field potential (LFP) recordings, with spatial resolution of 2-5 mm for ECoG and 0.1-1 mm for LFP, contrasting with the 20-30 mm of (high resolution) EEG.
We have shown that the model is readily simulated using standard finite difference schemes for PDEs, and validated one proposed scheme against an exact traveling front solution. We have performed parameter studies for the effects of system parameters on the speed of traveling fronts, and shown that correlations between somatic and dendritic coordinates in axodendritic connectivity patterns can strongly affect the speed of a wave. With an increase in the correlation parameter κ we see a slow down in the speed of the wave and the peak response shifts further away from the soma. For deployment in a neuroscience setting it is interesting to consider further work to optimize parameters of the model to fit real data. In this regard it is natural to extend previous work of Bojak and Liley, using particle swarm optimization, to generate EEG power spectra densities (PSDs) similar in shape to the ones observed in humans [33]. Such PSDs can be generated analytically using a noise-driven linear response theory as well as numerically thorough simulation of the full PDE model, after the inclusion of a noise source in Eq. (2). Importantly this would shed light on how real patterns of axodendritic connectivity can shape PSDs, as well as the known scaling of synaptic strength with distance from the soma that underlies so-called dendritic democracy [34]. Moreover, the model can also be used to generate a more direct measure  of EEG, by recognizing that the currents along the dendrites can be used to construct electromagnetic fields [9,35]. The transmembrane current I mem a at a point on the cable is equal to 1/r a ∂ 2 V a /∂x 2 , where r a is the (constant) specific resistance per unit length for currents flowing along the dendrite (and r a = 4R a /(π d 2 ), where d is the diameter of the dendrite and R a its axial resistivity). If the extracellular medium is assumed to be homogeneous, purely resistive, and infinite in extent with conductivity σ e , then the potential arising from the transmembrane currents is (in the quasisteady approximation to Maxwell's equations) given by where d (x, x , r, r ) = (x − x ) 2 + |r − r | 2 is a simple notion of distance. Indeed might give a more apt description of the fields observed using laminar implanted electrodes [36]. One final generalization of the work in this paper would be to treat more realistic models of dendrites. From a morphological perspective this could be incorporated by replacing the Green's function of the simple cable by one for a branched system. This can be constructed analytically using the sum-over-trips formalism, that itself can also cope with quasiactive membrane models (describing the resonance properties of dendrites associated with certain nonlinear ionic currents) [37]. All of the above are topics of ongoing research and will be reported upon elsewhere.

APPENDIX A: NUMERICAL METHOD
Here, we describe the numerical scheme for evolving a single population with only excitatory interactions. This allows us to suppress the indices a, b, although the extension to include inhibition is straightforward by generalizing the scheme below (and using the index notation of Sec. II). We note that expanding the differential operator A 2 in Eqs. (17) and (19) gives rise to a 'negative diffusion' term, that, if not handled sensibly, could give rise to numerical instabilities. To circumvent this potential issue we write A 2 ψ in the coupled form A 2 ψ = Ay, where y = Aψ. After discretizing in space the model in two spatial dimensions given by Eq. (19) takes the form of a system of coupled first order ordinary differential equations that we can write as where

t ), and (b)
x is the first order backwards difference operator acting on x, and (c) xx and (c) zz are the second order central difference operators in the x and z directions, respectively. Here, the variables (Y 1 , ψ,Y 2 , g, V ) are interpreted as vectors on the corresponding (uniform) spatial meshes.
A similar approach may be taken for the three-dimensional model given by Eq. (17) and yields essentially the same equations with finite-difference operators along the z direction replaced by the Kronecker sum of the second order central difference operators along the r 1 and r 2 directions which we denote by (c) rr , and the further replacement Here, V 0 = V (0, r, t ). The time integration was performed using an adaptive solver from the JULIA package-DifferentialEquations.jl [31,32].