Response functions for electrically coupled neuronal network: a method of local point matching and its applications

Neuronal networks connected by electrical synapses, also referred to as gap junctions, are present throughout the entire central nervous system. Many instances of gap-junctional coupling are formed between dendritic arbours of individual cells, and these dendro-dendritic gap junctions are known to play an important role in mediating various brain rhythms in both normal and pathological states. The dynamics of such neuronal networks modelled by passive or quasi-active (resonant) membranes can be described by the Green’s function which provides the fundamental input-output relationships of the entire network. One of the methods for calculating this response function is the so-called ‘sum-over-trips’ framework which enables the construction of the Green’s function for an arbitrary network as a convergent infinite series solution. Here we propose an alternative and computationally efficient approach for constructing the Green’s functions on dendro-dendritic gap junction-coupled neuronal networks which avoids any infinite terms in the solutions. Instead, the Green’s function is constructed from the solution of a system of linear algebraic equations. We apply this new method to a number of systems including a simple single cell model and two-cell neuronal networks. We also demonstrate that the application of this novel approach allows one to reduce a model with complex dendritic formations to an equivalent model with a much simpler morphological structure.


Introduction
Neuronal cells have a distinctive structure which differentiates them from any other cell types. The most extended parts of many neurons are dendrites, and their morphological complexity has fascinated scientists since the exemplary work of Ramón y Cajal [3]. Organised in a network, neurons receive and integrate thousands of neuronal inputs via both chemical and electrical synapses located primarily on dendrites. With the development of sharp micropipette electrodes, dynamic properties of dendritic membranes started to be revealed through intracellular recordings, and in the late 1950s experimental work was complemented with the pioneering theoretical work of Wilfrid Rall on the application of cable theory to dendritic modelling. Rall's significant contribution to the topic of dendritic function is nicely summarised in the book of Segev et al. [15]. Recent experimental and theoretical/computational studies at a single cell level reinforce the fact that dendritic morphology combined with membrane properties plays an important role in dendritic integration (two books, edited by Stuart et al. [16] and Cuntz et al. [7], give informative overviews from both an experimental and a theoretical perspective). An additional level of complexity associated with synaptic connectivity needs to be taken into consideration when dynamics of neuronal networks, rather than single cell dynamics, are investigated.
The dendritic membrane of various types of neurons is known to be equipped with voltage-gated ion channels, nonuniformly distributed throughout dendritic arbours and often demonstrating nonlinear dynamics. Many models of neuronal cells with retention of complex dendritic forma-tions are built by combining the linear (passive) properties of dendrites together with nonlinear (active) dynamics of ion channels. At the level of a single cell or at the network level, such models are restricted to being solved only by numerical methods, based on a compartmental approach [14]. Although the nonlinear properties of voltage-gated ion channels contribute considerably to neuronal input-output relations, it is important to recognise that the purely passive or resonant (quasi-active) properties of dendritic membranes provide the fundamental core for signal filtration and integration. Resonant dynamics of dendritic membrane are usually associated with the hyperpolarisation-activated I h current and, from a mathematical perspective, can be described by linearising channel kinetics [9][10][11].
Here we focus on a network of neuronal cells with purely passive or resonant membrane dynamics coupled by dendrodendritic electrical synapses, also known as gap junctions. Gap junctions are mechanical and electrically conductive links between adjacent neuronal cells that permit direct electrical connections between them. Having been first discovered at the giant motor synapses of the crayfish in the late 1950s, gap junctions are now known to be expressed in the majority of cell types in the brain [8;13]. Using the cable theory approach for modelling dendritic arbours, the response of an entire dendro-dendritic gap junction-coupled neuronal network to any injected current can be represented by a response function. This response function, often referred as a Green's function, describes the voltage dynamics along a network structure in response to a Dirac delta pulse applied at a given discrete location. One of the methods for constructing the Green's function, the so-called 'sum-over-trips' approach, is built on a path integral formulation and was originally proposed by Abbott et al. [1;2] for passive dendrites of a single cell and then generalised by Coombes et al. [6] for resonant membranes and Timofeeva et al. [19] for a neuronal network. This method calculates the response function as a convergent infinite series solution consisting of terms with various trips (paths) on a given branching structure and the associated coefficients obtained by the sum-overtrips rules. It has been shown at the single cell level that although convergence of this method is fast for simplified dendritic structures, the number of trips to guarantee a small convergence error for real morphologies might be large and have a strong effect on computational efficiency [4]. Here we propose an alternative method for calculating the Green's function on a neuronal network coupled by dendro-dendritic gap junctions. This new method, named as a method of local point matching, is inspired by the sum-over-trips approach and utilises the trip coefficients of that method, but avoids the construction of any trips. Instead, the new method is based on the construction of a linear system of algebraic equations and therefore leads to compact solutions without an infinite number of terms. In Sect. 2, we introduce the network model for gap junction-coupled neurons. Each neuron in the network consists of a soma and a dendritic arbour. Cellular membrane dynamics are modelled by a resonant electrical circuit. In Sect. 3, we develop a new method of local point matching from the generalised form of the sum-over-trips approach [19] for constructing the Green's function for an arbitrary network. Applications of this new method are demonstrated in Sect. 4. We start with a simple single cell model consisting of a soma and dendrite and then move to a two-cell simplified network and, finally, to a more complex tufted network. Not only do we apply the local point matching method for constructing the Green's functions in each case, but also use it to reduce the full two-cell tufted network model to an equivalent and much simpler model. The last two aforementioned sections include the key results and skip some mathematical details on the derivation of analytical results. We refer the interested reader to "Appendix" for detailed mathematical derivations. Finally, in Sect. 5, we provide a discussion of our results, as well as possible extensions of this work.

The model
We consider a network of neuronal cells. Each cell consists of an arbitrary structure of a dendritic morphology and a lumped soma, and cells in the network are connected by gap junctions (see an illustrative example for two cells in Fig. 1).
The transmembrane voltage V i (X, t) on an individual branch i of each cell is governed by the following set of equations: These equations provide a more general case of the linear cable theory with the cell membrane modelled by the socalled 'LRC' (or resonant) circuit instead of the 'RC' (or purely passive) circuit. The resonant circuit for each branch is described by the specific membrane capacitance C i , the resistance across a unit area of passive membrane R i and an inductance L i in series with a resistance r i . The presence of an inductive path in the circuit is the result of the linearisation of channel kinetics (in this case with a single nonlinear gating variable) responsible for subthreshold oscillatory behaviour around the steady state [6;9-11]. The constants D i and τ i in Eq. (1) can be found in terms of the electrical parameters of the cell membrane as where a i is the diameter and R a,i is the specific cytoplasmic resistivity of branch i. The term I inj,i (x, t) models an external current applied to this branch. The dendritic structure of each cell is attached to an equipotential soma of the diameter a S modelled by the 'LRC' circuit with the parameters C S = C soma πa 2 S , R S = R soma /(πa 2 S ), L S = L soma /(πa 2 S ) and r S = r soma /(πa 2 S ). Moreover, individual branches of different cells can be connected by gap junctions described by a conductance parameter g GJ .
Equations (1) and (2) for each dendritic segment must be accompanied by additional equations describing the dynamics of voltage at the two ends of a segment. If the proximal (X = 0) or distal (X = L i ) end of a branch is a branching node point, the continuity of the potential across a node and Kirchhoff's law of conservation of current are imposed. For example, boundary conditions for a node shown in Fig. 1 take the form: where r a, j = 4R a, j /(πa 2 j ) is the axial resistance of branch j. If a branch terminates at X = L i , we have either a no-flux (a closed-end) boundary condition or a zero value (an open-end) boundary condition A lumped soma can be treated as a special node point with the somatic membrane voltage V S (t) and the following set of equations which imposes special boundary conditions on the proximal ends of branches connected to the soma: where the sum in Eq. (8) is over all branches connected to the soma. If the branches of two cells are coupled by a gap junction, the location of this coupling can be treated as a special node point on an extended branching structure. This gapjunctional (GJ) node requires the following set of boundary conditions (given here with an assumption that it is placed at X = 0): where m − and m + (n − and n + ) are two segments of branch m (branch n) on the left and right from a gap junction (see Fig. 1), respectively. The expressions in (10) reflect continuity of the potential across individual branches m and n, and Eqs. (11) and (12) enforce conservation of current. The whole network in Fig. 1 can be viewed as a graph structure (which can be cyclic) with different types of nodes: a terminal, a regular branching node, a somatic node or the GJ node. The voltage dynamics along the network structure are described by linear equations, and therefore, the model's behaviour can be studied by constructing the network response function known as the Green's function, G i j (X, Y ; t). This function describes the voltage response at the location X on branch i in response to a Dirac delta pulse applied to the location Y on branch j at time t = 0.

Method of local point matching for finding the Green's functions
Our method is based on and developed from the sum-overtrips approach for calculating the Green's function on a network of electrically coupled neuronal cells [19]. Considering a network of cells as a single extended graph structure with labelled branches {1, 2, . . . , i, . . . , k, . . . , j, . . . }, the generalised sum-over-trips framework allows one to construct the Green's function for the whole structure in the  (2) with initial conditions V i (X, 0) = 0 and I i (X, 0) = 0, we obtain an ordinary differential equation for each branch i: where Considering an injected current in the form of a Dirac delta pulse and rescaling each branch k of the network by its own characteristic function γ k (ω) as L k = γ k (ω) L k , it is possible to derive (see [6;19]) that the Green's function on a scaled network (x = γ i (ω)X , y = γ j (ω)Y ) takes the form of an infinite series expansion where f (x) = e −|x| and L trip (x, y; ω) is the length of a trip along the network structure that starts at the point x = γ i (ω)X on branch i and ends at the point y = γ j (ω)Y on branch j. The trips coefficients A trip (ω) in (15) are chosen according to the following set of rules: • Somatic node: A trip (ω) is multiplied by a factor 2 p S,k (ω) or 2 p S,k (ω) − 1 (see Fig. 2b), where • GJ node: A trip (ω) is multiplied by a factor p GJ,n (ω), and R GJ = 1/g GJ . • Terminal: A trip (ω) is multiplied by +1 for the closed-end boundary or by −1 for the open-end boundary condition.
We refer the reader to [19] for a detailed summary of the generalised sum-over-trip method and the trip coefficients. Next, we provide a description of the main steps behind the derivation of the new method of local point matching together with the algorithmic summary of this method, the detailed derivation of which can be found in "Appendix 1". Note that ω is omitted for compactness from this point. All trips terminated at point y can be divided into two classes separated by the direction of the last part of the trip. Placing two points v j and w j on segment j as shown in Fig. 3, we consider one class which includes the trips with L trip (x, v →y j ) approaching y from the left (named as J v j ) and the other class which includes the trips with L trip (x, w →y j ) approaching y from the right (named as J w j ). Without constructing the actual trips, it is possible to show that all trips ending at y, named as J y and from (15) having the form can then be found as a linear combination of the unknown functions J v j and J w j belonging to these two classes. Likewise, we can partition trips on all other branches by placing a pair of points (v k , w k ) on each segment k and introducing two classes of trips J v k and J w k (see Fig. 4). Each unknown function J v k can then be written as a linear combination of This leads to a linear system of 2N algebraic equations for all unknown functions J v k and J w k defined on each segment k = 1, . . . , N , where N is the number of dendritic segments in the network. Solving this linear system for J v k and J w k , we can then find the unknown function J y and, as a result, the Green's function G i j (x, y; ω).

Summary of method
Here we summarise the main steps of an algorithm for constructing the compact Green's functions in the Laplace domain for an arbitrary neuronal network and refer the reader to "Appendix 1" for the detailed derivation of this method.
1. The physical length L k of each branch k is scaled by its own characteristic function γ k (ω) given by Eq. (14). 2. Place a pair of points (v k , w k ) on each segment k (see Fig. 4). Assume that v k and w k are placed infinitesimally close to both ends of the branch. Trips from v k and w k can move only towards each other (see red vectors in Fig. 4). erwise, an additional term a ik f (x) must be included in the linear combination, where a ik is a coefficient for a trip passing from segment i to segment k). The function J w k in Fig. 4 depends on a linear combination of the terms with J v k and J w n−1 . The constructed linear combinations for the unknown functions J v k and J w k include trip coefficients a nk for trips passing from segment n to segment k and trip coefficients a kk for trips reflecting at the end points of segment k. These coefficients are obtained from the sum-over-trips rules summarised in Fig. 2. 3. Solve the constructed linear system of algebraic equations and therefore find J v j and J w j for a pair of points (v j , w j ) placed on segment j which includes point y, and take the inverse Laplace transform (InvLT) of If point y is located at a node (i.e. y = 0 or y = L j ), due to the continuity of the potential at the boundaries the method can be easily applied by initially, assuming that y is placed on segment j slightly away from this node and, after the Green's function is constructed, considering that y = 0 or y = L j . A similar approach can be used if point x is also located at one of the nodes.
Note that spatially extended neurons coupled by gap junctions into an arbitrary neuronal network might develop a graph structure with cycles, and our method of local point matching (as well as the original sum-over-trips method) can support such structures.

A soma and dendrite model
Here we consider a simple model of a dendrite with a lumped soma attached to it at x = 0 (see Fig. 5). We assume that Somatic voltage profile in response to a stimulus I chirp (t) with parameters ω chirp = 0.003, A chirp = 0.2 nA a length of the dendrite, L, is scaled by its characteristic function γ (ω), i.e. L = γ (ω) L . If the dendrite is terminated with a closed-end boundary condition (i.e. has a factor +1 at the terminal), a system of linear equations for J v and J w corresponding to a pair of points (v, w) takes the following form where p S can be found from (17) and (18) as Solving the system, we can find that and then obtain J y as and finally, the Green's function in the Laplace domain This compact solution for the Green's function is equivalent to a solution in the form of an infinite series expansion obtained by using the sum-over-trips method [18]: If the output is measured at the soma (x = 0), the compact Green's function takes the form and in the case of the somatic stimulation, it is simply .
(30) Figure 6a shows a profile of the somatic Green's function given by (30). In Fig. 6b, we plot a somatic voltage profile in response to a chirp stimulus The model can be easily modified for the case of a semiinfinite dendrite. Assuming L → ∞, we have f (L) → 0 and from Eqs. (24) and (25) Then J y in (27) takes the form Other parameters as in Fig. 6 and the somatic Green's function simply becomes Resonant dynamics of the model can be characterised by a preferred frequency Ω * at which the Green's function has its maximum. Figure 6b clearly shows resonant behaviour of the system maximising the voltage response for particular frequencies. In Fig. 7, we plot a preferred frequency as a function of a dendritic length L when x and y are placed at the soma. This dependence is obtained as a solution of the implicit equation The plot demonstrates a nonmonotonic trend with a minimal value within a realistic range of dendritic lengths indicating a considerable effect of dendritic extents on the model's dynamics.

A two-cell simplified network
Here we demonstrate how our method can be applied to a twocell network of either identical or nonidentical cells coupled by a dendro-dendritic gap junction. In each case, we obtain the compact solutions for the Green's functions, Eqs. (40)-(43) for the two-cell identical network and Eqs. (49)-(56) for the two-cell nonidentical network, which can inform us about the roles of individual parameters on the network dynamics.
We start by considering a model of two identical cells, either of which consists of a soma and N attached semiinfinite dendrites as shown in Fig. 8. We assume that the biophysical properties of all dendritic segments are the same and that the physical lengths are scaled by the characteristic function γ (ω) given by (14). The gap junction is located at some distance L GJ away from the cell bodies. We assume that this network can receive stimuli in four different locations mimicking distal (y 1 and y 2 ) and proximal (y 3 and y 4 ) inputs.
Points of output x 1 (for Cell 1) and x 2 (for Cell 2) are located between either soma and the gap junction. Using our method, we can construct a linear system of algebraic equations for the functions J a , J b , J y and J w in the case of placing output at x 2 (see Fig. 8): This system of equations can be easily solved analytically (see "Appendix 2"). The Green's functions for four individual inputs for Cell 2 can then be found as and Here p S and p GJ can be found from Eqs. (17) and (19), respectively. As the cells are identical and due to the symmetry of the input locations, the corresponding Green's functions for Cell 1 can be easily obtained from Eqs. (40)-(43). In Fig. 9, we plot the Green's functions at the soma of each cell (x 1 = 0 and x 2 = 0) in response to a stimulus y 3 = 0 applied to Cell 1. Note that Eqs. (40) and (41) are equivalent to the solutions for the Green's functions in the form of an infinite series expansion found using the 'method of words' in [19]. Truncated series solutions with the index of summation n ≥ 10 match the compact Green's functions obtained here (not shown).
Assuming that two distal inputs y 1 and y 2 are applied at equal distances from each soma (y 1 = y 2 > L GJ ), the Green's function for each soma is identical: Similarly, for the case of two proximal inputs y 3 and y 4 placed at the same distance away from each soma (y 3 = y 4 < L GJ ), the somatic Green's function for each cell has the same form: Both solutions are independent of g GJ and L GJ and have the form of Eq. (34) for the soma and dendrite model, i.e. without the presence of the gap junction. Using our method, we can also construct analytical solutions for the Green's functions for a network of two nonidentical cells. The response functions at the soma of Cell 2 take the forms and using symmetry, the response functions at the soma of Cell 1 can be found as  Fig. 9, except r 1 = 100 · cm 2 . The dendritic parameters of Cell 2: a 2 = 0.4 μm, C 2 = 1 μF · cm −2 , R 2 = 20,000 · cm 2 , R a,2 = 150 · cm, and r 2 → ∞ (i.e. passive dendritic membrane). Both somas are passive where, for k = 1, 2, is the characteristic function of the membrane of Cell k, L k is the distance between the gap junction and the soma of Cell k, and p GJ,k is given by Eq. (19). Using Eqs. (49)-(56), we can investigate how the strength and location of the gap junction affect the dynamics of the two-cell model. Here, we consider that a stimulus is applied to the soma of Cell 1 and construct a map for the preferred frequencies Ω * 1 and Ω * 2 at the soma of Cell 1 and Cell 2, respectively. This map is shown in Fig. 10. In this case, Cell 2 is assumed to be purely passive, and Cell 1 has a passive soma with resonant dendrites. The map indicates that the location of a gap junction plays a significant role in the dynamics of the network, unless the coupling is weak. Moreover, the initially passive soma of Cell 2 starts to demonstrate a resonant behaviour imposed by Cell 1 even for weak coupling.
Often it is difficult to measure experimentally locations and strengths of gap junctions in real neuronal networks. Knowledge of the inverse map from a pair of preferred frequencies (obtained from somatic sub-threshold stimulations) to (L GJ , g GJ ) might provide estimates for gap-junctional parameters. However, the map Ψ is neither surjective nor injective (see, for example, Fig. 11 for a network of two resonant cells showing that the system may demonstrate the same resonant behaviour for two different gap-junctional locations, proximal and distal, and identical coupling strengths) making it mathematically impractical to obtain Ψ −1 . At the same time, if a constraint on locations of gap junctions is imposed (e.g. proximal or distal), this may lead to a one-to-one correspondence between (L GJ , g GJ ) and (Ω * 1 , Ω * 2 ) and therefore assists in the estimation of

A two-cell tufted network
Now we consider a more realistic neuronal network consisting of two identical tufted or mitral cells. Each neuron has a soma attached to N dendritic branches, one of which is the primary dendrite with the tuft spanning from its end. Two cells are coupled in their tufts by dendro-dendritic gap junctions (see Fig. 12a). As in the previous model, we assume that the biophysical properties of all dendritic segments are the same and that the physical lengths are scaled by the characteristic function γ (ω). We consider that each cell has n T segments in its tuft, and n GJ of them possess identical single gap-junctional points located l 0 away from the end of the primary dendrite. The primary dendrite of each cell has the length L, while the other branches are semi-infinite. For simplicity, we consider that the membrane of both cells is purely passive (i.e. γ 2 (ω) = (τ −1 + ω)/D); however, it is straightforward to generalise it for the resonant case.
Although it is possible to use our method and construct the compact Green's functions for this tufted network, we will first demonstrate that there exists an equivalent reduced model with the simplified structure shown in Fig. 12b for which the compact solutions will then be constructed. Here we consider a reduction in the full model when external inputs cannot apply to any of the tufted dendrites. Notating the reduced network with the symbol we constrain the reduced model to have where p S constitutes part of the definition of a trip coefficient for the somatic node (see Fig. 2b), and p D is a branch factor of the primary dendrite defined as in (16) and constitutes part of the definition of a trip coefficient for the branching node (see Fig. 2a). Equations (63) and (64) force the length of the primary dendrite and the location of a single gap junction in the reduced model to be the same as in the full tufted model. Placing y 0 on the primary dendrite of Cell 1, both models are equivalent if Using our method of local point matching, we can write down a system of algebraic equations for the full tufted model: It is possible to prove that Eq. (67) holds and the systems in Fig. 12a, b are equivalent when, in addition to constraints (63)-(66), p T = n GJ p T , R GJ = R GJ /n GJ and z = n GJ z, giving which constitutes part of the definition of a trip coefficient for the GJ node (see Fig. 2c). A detailed proof of model reduction is given in "Appendix 3". Our method for constructing the compact Green's functions can then be simply applied to the reduced model shown in Fig. 12b. A detailed derivation of the solutions for the model with a stimulation applied at y 0 in Cell 1 and the output points x 0 placed in both cells is given in "Appendix 4". Assuming that y 0 is placed at the soma of Cell 1, the Green's functions at each soma have the following forms where For investigating the effect of gap junctions from the tufted regions of the cells on the model's behaviour, we define a coupling ratio (CR) as Using Eqs. (87) and (88), we compute and plot in Fig. 13 a map of CR for various values of conductance g GJ and location l 0 of the gap junctions in the tuft. This map can be compared with the CR map obtained earlier in [12] for two mitral cells coupled by distal gap junctions. Note that the map in [12] is obtained by brute-force numerical simulations of a computational model with a similar, but not identical, structure to our two-cell model. Using our method of local point matching, we can also prove that there exists an equivalent reduced model for the full tufted model with external inputs applied to the tufts (instead of the primary dendrites). We consider that any tuft dendrite k can receive a Dirac delta pulse at the location y k away from the branch point with the primary dendrite. This tuft dendrite can be either with or without a gap junction. In the equivalent reduced model shown in Fig. 14, we consider two possible inputs corresponding to the location of y k , namely the input y 1 applied to the branch without a gap junction and y 2 applied to the branch with a gap junction. It is possible to show (see "Appendix 5" for details) that the Green's function of the full tufted model for a given input y k can be found knowing the Green's function for the equivalent reduced model as for the input y k applied to the branch without a gap junction, and for the input y k applied to the branch with a gap junction.
Here, the reduced model is constructed in such a way that the stimuli in the full and reduced models are located at the same distance away from the primary dendrites, i.e. y 1 = y k and y 2 = y k . The point x 0 (0 ≤ x 0 ≤ L) is located on the primary dendrite of either of two cells. The reduced model shown in Fig. 14 can be constructed from the full tufted model using a number of constraints specified in (164)-(168) and (174). The Green's functions G for each cell can then be found by our method of local point matching and used in (90) and (91) for finding the Green's functions for the full tufted model. In the case of multiple inputs applied to the tuft dendrites, the Green's function for each cell can be found by summing individual Green's functions for each input. Assuming that all tuft dendrites of both cells receive identical inputs located at a distance y away from the primary dendrite of each cell, we obtain in this special case

Discussion
In this paper we have presented a novel method for calculating the Green's functions for arbitrary neuronal networks with a passive or resonant cell membrane coupled by dendrodendritic gap junctions. This method provides an alternative and complementary approach to the generalised sum-overtrip method [19]. Importantly, our new approach avoids the construction of an infinite number of trips and, being based on the construction of a linear system of algebraic equations, provides exact expressions for the network response function in the Laplace (frequency) domain without any issues of computational convergence. We have applied this new method of local point matching to a simple single cell model and two-cell neuronal networks (simplified and with tuft dendrites). Its application to the tufted network has also allowed us to reduce it to an equivalent network, but with a much simpler morphological structure. We have also illustrated that knowledge of the exact compact expressions for the Green's function can provide important insights into the role of individual variations in cell parameters on the model's dynamics.
There are a number of natural extensions of the work in this paper. One is an application to more realistic network geometries with more than just two cells, given that a computational implementation of the method of local point matching can provide a fast realisation of the Green's function for the whole network. Having a complex network of multiple cells with a graph structure consisting of N dendritic segments, we need to construct and solve a linear system of 2N equations only once to find all unknown J v k and J w k functions. We can then simply construct the functions J y for each dendritic segment to obtain G i j (X, Y ; ω). Note that the point X can be placed on each dendritic segment before constructing a system of linear equations for J v k and J w k . Switching off all X points except one on branch i in the solution for J y allows one to find the Green's function for the entire network straight away. The numerical inverse Laplace transform to obtain G i j (X, Y ; t) is the only procedure in which a computational approximation appears. As has been previously pointed in Sect. 4.2, knowledge of a map from the preferred frequencies to the system's parameters for a reconstructed neuronal network combined with subthreshold electrophysiological data might provide some estimates for important network's parameters and additional work is required in this direction. Another possible extension is to incorporate active properties in dendrites and somas of cells in a network and analyse the propagation of dendritic action potentials as well as somatic spiking dynamics. The spike-diffuse-spike (SDS)-type model [5;17] can be utilised for that, as although the voltage-gated channels in the SDS framework are modelled by piecewise linear instead of nonlinear dynamics, it has been shown that the speed of wave propagation in the SDS model is in excellent agreement with a more biophysically realistic nonlinear model [20]. Both these extensions will be reported on elsewhere.
Note that ω is omitted for compactness from this point. The trips in (95) are divided into two groups: the trips that are passing through v j just before reaching y and the trips that are passing through w j just before reaching y (see Fig. 3). L trip (x, v →y j ) introduced in (96) defines a length of a trip which moves in the direction of y and ends at v j before reaching y. Similarly, L trip (x, w →y j ) defines a length of a trip which moves in the direction of y and ends at w j before reaching y (shown in red in Fig. 3). A v j and A w j are the trip coefficients corresponding to the trips to v j and w j and obtained using the rules of the sum-over-trips method.
As v j is placed close to one end of the segment, we have L trip (x, v j ) = L trip (x, v →y j ), and therefore, we introduce Similarly, as w j is located close to the other end of the segment, we introduce Then simplifying the notations as J i j (x, y) = J y , J i j (x, v j ) = J v j and J i j (x, w j ) = J w j we can rewrite Eq. (96) in the form As both points v j and w j are placed infinitesimally close to individual ends of segment j of length L j , we consider that v j = 0 and w j = L j , and therefore, Eq. (99) can be rewritten as If the point y is located on a semi-infinite branch and w j is placed on the side towards infinity, then |w j − y| → ∞ giving f (w j − y)J w j = 0. Following similar steps, if placing two points v k and w k on each segment k infinitesimally close to each end, we can define functions J v k and J w k which can be written in terms of functions J v n and J w n associated with points v n and w n from all branches connected to a single node. For example, given a node with K segments and K pairs of points (v k , w k ) (see Fig. 4), the function J v k can be found as where a nk is a coefficient for a trip passing from segment n to segment k, a kk is a coefficient for a trip reflecting at one of the ends of segment k and L n is the scaled length of branch n. Equation (101) can be constructed for any node branches of which do not include point x. If x is located on branch i connected to a node in consideration (0 < x < L i ), an additional term representing a trip from x to v k needs to be added: For a given network, the number of functions J v and J w is equal to the degree sum of the corresponding graph and a system of linear equations for unknown J v and J w can be written using Eqs. (101) and (102). It is possible to show that this system of equations is linearly independent and therefore has a unique solution. Solving it, we can find J v j and J w j and obtain J y = J i j (x, y) from Eq. (100). The Green's function G i j (x, y) can then be calculated from Eq. (93) as If both x and y are located on the same branch, equation for J y takes the following form (instead of Eq. 100):
Using Eqs. (38), (39) and (106), we can obtain that where function F is defined in (46). Placing the stimulus to the point y 1 (see Fig. 8), we have with and the solution G 2 (x 2 , y 1 ) can be found as Similarly for y 2 , we can show that Solutions G 1 (x 1 , y 2 ) and G 1 (x 1 , y 1 ) can be easily obtained using the symmetry of the system. Note that for the case of two symmetrical inputs located after the gap junction (i.e. y 1 = y 2 = y), we have When the stimulus is applied to the location y 3 , we can get that and therefore obtain To find G 2 (x 2 , y 4 ), we notice that if two symmetrical inputs (y 3 = y 4 = y) are applied to the system, then using Eqs.

Appendix 3: Model reduction in a two-cell tufted network
In the case of the model with the input/output points placed in the primary dendrite, it is straightforward to demonstrate using our algorithm that n T − n GJ branches having no gap junctions can be easily merged into a single branch. From Eqs. (75), (76), (84) and (85) we obtain that Then due to Eq. (67), we get and then J y = J y .
Comparing Eqs. (74) and (83) Placing y 0 and x 0 on the primary dendrites of Cell 1 and Cell 2, respectively, we can write that As all dendrites in the tuft are identical, Eq. (135) is simply p T = n GJ p T . Using this relationship in (149) and having x 0 = y 0 = 0, we obtain the somatic Green's function for Cell 2 (Eq. 87).
As the cells are identical and the system is symmetrical, the Green's function for Cell 1 (where the input is placed) is equal to the Green's function for Cell 2 if both x 0 and y 0 are placed on Cell 2. For this case J y 0 takes the form Then G 2 (0, 0) + G 1 (0, 0) = 1 2Dγ [J w f (L) + J y f (0) which gives Eq. (88). Note that by setting p D = 1/2 which gives δ = 0, p T = 1/2, θ = f (L), we can recover a case of two simplified identical cells (without the tuft). Equation (87) gives and this is exactly Eq. (42) with x 2 = 0, y 3 = 0 and L GJ = L + l 0 .