Skip to main content
Advertisement
  • Loading metrics

Ephaptic coupling in white matter fibre bundles modulates axonal transmission delays

  • Helmut Schmidt ,

    Roles Conceptualization, Formal analysis, Investigation, Methodology, Visualization, Writing – original draft, Writing – review & editing

    hschmidt@cbs.mpg.de

    Affiliation Max Planck Institute for Human Cognitive and Brain Sciences, Leipzig, Germany

  • Gerald Hahn,

    Roles Conceptualization, Writing – review & editing

    Affiliation Center for Brain and Cognition, Computational Neuroscience Group, Department of Information and Communication Technologies, Universitat Pompeu Fabra, Barcelona, Spain

  • Gustavo Deco,

    Roles Conceptualization, Funding acquisition, Writing – review & editing

    Affiliations Center for Brain and Cognition, Computational Neuroscience Group, Department of Information and Communication Technologies, Universitat Pompeu Fabra, Barcelona, Spain, Institució Catalana de la Recerca i Estudis Avançats (ICREA), Barcelona, Spain, School of Psychological Sciences, Monash University, Melbourne, Australia

  • Thomas R. Knösche

    Roles Conceptualization, Funding acquisition, Methodology, Supervision, Writing – review & editing

    Affiliations Max Planck Institute for Human Cognitive and Brain Sciences, Leipzig, Germany, Technische Universität Ilmenau, Institute of Biomedical Engineering and Informatics, Ilmenau, Germany

Abstract

Axonal connections are widely regarded as faithful transmitters of neuronal signals with fixed delays. The reasoning behind this is that extracellular potentials caused by spikes travelling along axons are too small to have an effect on other axons. Here we devise a computational framework that allows us to study the effect of extracellular potentials generated by spike volleys in axonal fibre bundles on axonal transmission delays. We demonstrate that, although the extracellular potentials generated by single spikes are of the order of microvolts, the collective extracellular potential generated by spike volleys can reach several millivolts. As a consequence, the resulting depolarisation of the axonal membranes increases the velocity of spikes, and therefore reduces axonal delays between brain areas. Driving a neural mass model with such spike volleys, we further demonstrate that only ephaptic coupling can explain the reduction of stimulus latencies with increased stimulus intensities, as observed in many psychological experiments.

Author summary

Axonal fibre bundles that connect distant cortical areas contain millions of densely packed axons. When synchronous spike volleys travel through such fibre bundles, the extracellular potential within the bundles is perturbed. We use computer simulations to examine the magnitude and shape of this perturbation, and demonstrate that it is sufficiently strong to affect axonal transmission speeds. Since most spikes within a spike volley are positioned in an area where the extracellular potential is negative (relative to a distant reference), the resulting depolarisation of the axonal membranes accelerates the spike volley on average. This finding is in contrast to previous studies of ephaptic coupling effects between axons, where ephaptic coupling was found to slow down spike propagation. Our finding has consequences for information transmission and synchronisation between cortical areas.

Introduction

Signal processing and transmission in neuronal systems involves currents flowing across neuronal cell membranes. Due to the resistance of the extracellular medium, such transmembrane currents generate extracellular potentials (EPs), also called local field potentials (LFPs). The sources of EPs are synaptic currents, action potentials, calcium spikes and voltage-dependent intrinsic currents [1]. Neurons can therefore interact with their neighbours by changing the electric potential of the extracellular medium (and hence the membrane potential of their neighbours) without forming synapses. Such interaction is termed ephaptic interaction or ephaptic coupling [24]. Since EPs generated in the cortex are generally of the order of 100μV [5] and therefore small in comparison to neuronal threshold potentials, the influence of EPs on neural computation is often regarded as negligible. EPs can be measured with intracranial electrodes and are used as a proxy for the underlying neuronal activity [69].

Seminal experiments by Katz and Schmitt [10], Rosenblueth [11], Arvanitaki [2] and Marrazzi and Lorente de Nó [12] have demonstrated that action potentials travelling along parallel axons can interact with each other if the extracellular medium is highly resistive. They demonstrated that action potentials with an initial offset would resynchronise, and also slow each other down. Furthermore, action potentials could be initiated in passive axons by action potentials travelling in a nearby axon. Several studies have reproduced these effects using computational models [1320]. However, the experimental setup is such that the axons are placed into a highly resistive medium (either paraffin oil [10], or moist air [11]) in comparison to the intracellular medium, and the computational models assume that axons are embedded within a finite-sized extracellular medium. The latter would be justified by the presence of epineuria or perineuria, which is tissue restricting the extracellular space around axons in the peripheral nervous system. Both, however, are unlikely scenarios for axonal fibre bundles within the brain: the extracellular medium is only about three times more resistive than the intracellular medium, and axons in the central nervous system are not wrapped by epineuria and perineuria that would justify the ‘cables within a cable’ approach. For these reasons, the amplitude of extracellular potentials around spike carrying axons should be small, and ephaptic coupling should not play a significant role between individual pairs of axons within axonal fibre bundles in the brain. However, we hypothesise that collective interaction between multiple axons affects axonal signal transmission.

We test our hypothesis by introducing a modelling framework in which EPs modulate spike thresholds, and hence spike propagation velocities. We first determine the EPs generated by action potentials in single axons, which can be computed using the axial profile of an action potential (Fig 1A). The importance of computing the EPs lies in the fact that they perturb the membrane potential of a passive fibre (Fig 1B). This is then followed by the computation of EPs generated by spike volleys in fibre bundles (Fig 1C and 1D). As axon bundles contain millions of axons, we compute the cumulative effect of spike volleys at the macroscopic scale in axon bundles with diameters of several millimeters. The results of this analysis are used to build a spike propagation model, in which each spike travels with a velocity that is determined by structural parameters of the axon and the extracellular potential. This model is then coupled into a neural mass model to investigate in-silico the latency-intensity relationship of sensory stimuli, and the role of ephaptic coupling.

thumbnail
Fig 1. Computing the extracellular potential (EP) generated by a volley of spikes.

A: An action potential, as expressed by the membrane potential Vm along the axial dimension z, generates an EP that varies with z and the distance from the axon d. B: An action potential in an active axon perturbs the membrane potential of a passive axon via the EP. C: We consider spike volleys travelling along axonal fibre bundles, and D: infer from the EP the cumulative effect on the membrane potential of a passive axon.

https://doi.org/10.1371/journal.pcbi.1007858.g001

The spike propagation model that we propose in this article constitutes a strong simplification of biophysical models for spike propagation in (myelinated) axons. Rather than numerically computing the membrane potential along the axon, we asign each spike a position on an axon that changes in time according to its propagation velocity. In the absence of any external perturbations, the propagation velocity is constant along the axon and scales with its structural parameters. For example, the propagation velocity increases linearly with the axon diameter. The effect of extracellular potentials is then modelled by a coupling function that rescales the propagation velocity as a function of the EP at the spike position. The EP in the fibre bundle is computed based on core conductor theory. Each spike that travels along the fibre bundle is asigned a characteristic spatial profile, the length of which scales linearly with the propagation velocity. As this model is a strong simplification of standard biophysical models, we use a computationally feasible scenario where a synchronous spike volley interacts with itself as a test bed to calibrate the coupling term in the simplified model.

There is no direct evidence how spike propagation in axon bundles is affected by ephaptic coupling effects. We therefore present indirect evidence based on psychological experiments that investigate the relationship between the intensity of a sensory stimulus and the latency between stimulus presentation and the maximum response of the evoked potential. Such experiments have been performed for a range of sensory stimuli, including visual [21], auditory [2225], and nociceptive stimuli [26, 27]. For auditory stimuli, the first maximum (P1) is observed about fifty milliseconds after stimulus presentation, and the drop between low-intensity and high-intensity stimuli is of the order of ten milliseconds [23]. All neuronal signals, including sensory stimuli, have to pass through axon bundles to reach cortical areas. Consequently, we set up a model system in which an axon bundle is coupled to the Jansen-Rit neural mass model, which is capable to produce evoked potentials. A crucial assumption we make here is that sensory stimuli are coded as highly synchronised spike volleys, or sequences thereof. Such spike volleys could be generated by sensory neurons that encode rapid changes in sensory stimuli, or by the interaction between excitatory and inhibitory neurons. Computational studies have demonstrated that cortical circuits are capable of generating highly synchronised spike volleys with millisecond duration [28, 29].

Results

Extracellular potentials around single axons

First we computed the EPs generated by action potentials in single axons. We used the line approximation [3032], given that the diameters of axons are several orders of magnitude smaller than the diameter (or the general lateral dimensions) of axonal fibre bundles: (1) Here, ϕ is the EP, z is the axial dimension, d the distance from the axon, t is time, σi and σe are the intracellular and extracellular conductivities, a is the axon radius and V″(z) is the second derivative (curvature) of the membrane potential V. The EP is computed for different approximations of the spatial profile of an action potential, which include a piecewise linear and a piecewise quadratic approximation of spike profiles, but also spike profiles generated by a biophysical model [33] (Fig 2A–2C). The advantage of the piecewise approximation of the action potential profile is that the EP can be computed analytically (see the Methods section for details). The EP obtained from the biophysical model is computed numerically. For all the profiles we find that the maximum amplitude of the EP is of the order of microvolts (Fig 2D–2F), and at long distances d the EP decays with d−3 (Fig 2G–2I), akin to electric potentials of quadrupoles.

thumbnail
Fig 2. Spatial profiles of action potentials and their EPs.

Shown are A: the piecewise linear profile, B: the piecewise quadratic profile, and C: the profile of an action potential generated with the biophysical model. D-F: EPs corresponding to action potential profiles in A-C. G-I: Log-log plots of the EPs (absolute values) at z = 0. Black lines indicate decay with d−3. (The notch at d ≈ 0.3mm is due to a change of sign).

https://doi.org/10.1371/journal.pcbi.1007858.g002

Extracellular potentials in fibre bundles

To compute the effect of multiple action potentials in a fibre bundle, we assumed that a completely synchronous spike volley travels through the fibre bundle. The fibre bundle was arranged as a set of concentric rings of axons, as shown in Fig 3A. The reference point to compute the EP was set at the centre of the axon bundle. We computed the EP for an increasing number of spikes, beginning with six spikes in the innermost ring of axons, then 18 spikes in the two innermost rings, and successively increasing the number of rings in which all axons carry action potentials (Fig 3A). The maximum number of rings considered in this setup was 104, which corresponds to a fibre bundle diameter of 10mm if the diameter of the uniform axons is 0.5μm. This fibre bundle contains approximately 3 × 108 axons, similar to the number of axons in the corpus callosum [34].

thumbnail
Fig 3. EP at the centre of a circular axon bundle due to concentric spike volleys.

A: Microscopic cross-section of a fibre bundle, with spike-carrying axons marked in blue. B: Macroscopic extension of (a), with the active area (i.e. where axons carry spikes) marked in blue. C: Waveform of a spike (top), and the resulting spatial (axial) profile of the EP at the centre of the fibre bundle. D: Cross-sections of C.

https://doi.org/10.1371/journal.pcbi.1007858.g003

Increasing the active area (see Fig 3B for a macroscopic representation) yielded a longitudinal profile of the EP that saturated at large diameters (Fig 3C). Interestingly, the profile is approximately proportional to −V(z), with V(z) being the spatial profile of the action potential (Fig 3D). In the Methods section we demonstrate that this profile can be computed analytically, to a very good approximation, by the following expression: (2) Here, g represents the ratio between the axon diameter and fibre diameter (axon plus myelin), commonly referred to as g-ratio. The relative size of the volume occupied by fibres, the fibre volume fraction, is represented by the quantity ρ, and P is the radius of the fibre bundle.

Next, we investigated how the EP changes with the position of the reference point, i.e. the point in the cross-sectional plane at which the EP is computed (Fig 4A). We found that the amplitude and longitudinal profile remained nearly unchanged, even if the reference point is close to the surface, as shown in Fig 4B and 4C. More specifically, the decrease of the amplitude is less than ten percent when the reference point is moved from the centre of the fibre bundle to 0.8 bundle radii away from the centre. Closer to the surface, the drop in amplitude is more marked. Outside of the bundle, while moving the reference point further away from the centre the EP drops rapidly, and at sufficiently large distances the drop in amplitude is proportional to d−3. We take this as evidence that the EP at the centre of the bundle is characteristic for the EP across the entire cross-section of the fibre bundle, i.e. we assume the EP is uniform in the radial direction.

thumbnail
Fig 4. EP in fibre bundle with synchronous spike volley, subject to position of reference point.

A: The reference point is moved from the centre of the fibre bundle to a position outside of the fibre bundle. B: Waveform of a spike (top), and the resulting EP plotted against the longitudinal coordinate z and the distance of the reference point from the centre. C: Cross-sections of B.

https://doi.org/10.1371/journal.pcbi.1007858.g004

We consider spike volleys that engage all axons in the fibre bundle, which leads to EPs with amplitudes of order 100mV, as can be seen in Figs 3C and 4B. This is certainly an unphysiological scenario, since it is unlikely that all axons in a fibre bundle carry perfectly synchronised action potentials, and because such large EPs would certainly disrupt signal transmission in the participating axons. However, it is plausible that a (sufficiently synchronous) spike volley might engage one percent of the axons in the fibre bundle, in which case the amplitude of the EP would be of the order of 1mV.

Alternatively, one may consider a spike volley that is not perfectly synchronised, i.e. the spikes are distributed in space due to varying emission times. To illustrate the effect of such a spatial distribution, we draw spike positions randomly from a uniform distribution of varying width Δz. This spatial distribution can be associated with a temporal distribution via the relation Δz = vΔt, where v is the (intrinsic) propagation velocity of the uniform axons. In Fig 5 we show how increasing the active area affects the EP for different Δz. It can be seen that the maximum amplitude decreases with increasing Δz, and for sufficiently wide spike volleys the largest amplitude of the EP occurs near the edges of the spike volley instead of its centre.

thumbnail
Fig 5. Increasing the length of a spike volley attenuates the amplitude of an EP.

The EP is shown for varying bundle diameters and z. We steadily increase the width Δz of the spike volley from A: Δz = 0mm, to F: Δz = 50mm.

https://doi.org/10.1371/journal.pcbi.1007858.g005

A model for spike propagation

In addition to studying EPs generated by spike volleys in axonal fibre bundles, we are interested in the effect that EPs have on axonal signal transmission. Since the membrane potential is measured as the difference between intracellular and extracellular potential, a change of the EP implies a change of the membrane potential. For example, if the EP decreases, then the membrane potential increases, i.e. the membrane is depolarised. We assume that the EPs are not compensated by transmembrane currents, or if such processes occur, that these processes are too slow to be relevant for short spike volleys.

We begin the modelling procedure by setting up a fibre bundle with N axons, each of which has a diameter drawn from a shifted alpha distribution that was chosen to closely fit the results by Liewald et al. [35] (Fig 6A). For numerical purposes we set the number of axons N between 103 and 104. A realistic fibre bundle contains many more axons, likely by several orders of magnitude. Conceptually, each of our model axons therefore represents a large number of axons with identical properties, but evenly distributed across the cross-sectional area. The fibre bundle is also endowed with macroscopic properties, namely the length and radius of the fibre bundle.

thumbnail
Fig 6. Illustration of properties of the computational model.

A: Distribution of axon diameters sampled from a shifted alpha distribution to match experimental data [35]. B: Rastergram of spike volley generated at proximal end of fibre bundle. C: Rastergram of spike volley reaching the distal end of the fibre bundle. D: Distribution of delay times. E: Snapshot of the longitudinal profile of EP generated by a spike volley. F: The EP modulates the spiking threshold (Vthr) and therefore the delay Δt of action potential generation between two reference points (e.g. two consecutive nodes of Ranvier).

https://doi.org/10.1371/journal.pcbi.1007858.g006

To test the transmission properties of a fibre bundle, we set up a spike volley with spike times drawn from a uniform distribution. The spike times define when the action potentials are generated at the proximal end of the bundle (Fig 6B). The propagation of spikes along the axon is determined by a spike propagation model that is described in the next paragraph. The spike volley then reaches the distal end of the fibre bundle (Fig 6C). Due to the distribution of axon diameters, this process results in a distribution of transmission delays (Fig 6D). If the position of a spike is known, one can determine the EP generated by this spike. Since each model axon represents a large number of biological axons, we do not use the expression for single axons (Eq (1), but the one for the cumulative EPs generated by spike volleys (Eq (2)). The EP generated by a spike is thus the EP shown in Fig 3C, divided by the volume fraction occupied by the model axon. In this way, one can compute the spatial profile of the EP generated by a spike volley, see Fig 6E.

The spike propagation model tracks the position of an action potential along the fibre bundle, which is determined by the leading edge (rising phase) of the action potential. For the linear and quadratic approximations of the spike profile, the position is defined by the point where the membrane potential first deviates from resting potential. In the absence of perturbations by non-zero EPs, the velocity is constant along the fibre bundle. Therefore, the position of a spike can be tracked by multiplying the intrinsic velocity (determined by structural parameters of the axon) with the time elapsed since the spike was generated. The velocity of a spike is also determined by a putative spike threshold, which might be interpreted within a spike-diffuse-spike framework [19]. It has been demonstrated, using some simplifying assumptions, that the spike threshold Vthr can be related to the activation delay Δt between two subsequent nodes of Ranvier by some nonlinear function, and therefore to the velocity of a spike [19]. In the presence of EPs, the membrane potential, and therefore the propagation velocity, is perturbed. The perturbation of the membrane potential can be intepreted as a perturbation of the spike threshold. If the membrane is depolarised (hyperpolarised) by the EP, then the spike threshold is effectively lowered (raised). For simplicity, we assume a linear relationship between Vthr and Δt (Fig 6F). This results in the following relationship between the perturbed propagation velocity v and the EP: (3) Adjusting the prefactor γ allows for the calibration of the spike propagation model with a more detailed biophysical model, which we demonstrate next.

Model calibration with biophysical model

An axon can be regarded as a core-conductor, and the spatio-temporal evolution of the membrane potential V can be described by the following cable equation: (4) The term on the left hand side describes capacitive trans-membrane currents. The first term on the right hand side describes changes in axial currents inside the axon, and the second and third term on the right hand side describe resistive currents across the axonal membrane and the myelin sheath. The second term represents passive currents, and I(V) represents voltage-gated currents, as described by the Hodgkin-Huxley model. Eq (4) describes the scenario when the extracellular potential ϕe is zero, in which case the membrane potential V equates the intracellular potential ϕi. If the EP is not zero, then the appropriate equation for the resulting membrane potential V = ϕiϕe is (5) The EP now affects the resistive currents (second and third term on the r.h.s.), as well as the capacitive current (fourth term on the r.h.s.).

We focus the calibration effort on a computationally feasible example. We consider a synchronous spike volley in a bundle composed of identical axons. The spike volley consists of spikes in one percent of all fibres. As all spikes are identical, Eq (36) is representative for all spike-carrying axons. The EP is calculated numerically at each time step from the resulting profile of the membrane potential using Eq (2). We vary the axon bundle diameter and record the change in the propagation velocity for the biophysical model and the spike propagation model (Fig 7).

thumbnail
Fig 7. Comparison of the spike propagation model with a biophysical model.

A synchronous spike volley slows down as a result of ephaptic coupling in fibre bundles with identical axons. The relative change of the propagation velocity varies with the bundle diameter.

https://doi.org/10.1371/journal.pcbi.1007858.g007

To fit the spike propagation model to the biophysical model, we adjust the coupling parameter γ and the standardised spike profile. Specifically, we adjust the length of the spike profile and the position of the maximum. The parameter γ determines the amount of relative slowing down that is observed in both models, whereas the profile of the model spikes determines how this effect changes with varying bundle diameters (Fig 7).

Effect of extracellular potentials on transmission delays

The spike propagation model allows us to test the consequences of ephaptic coupling via EPs in a macroscopic fibre bundle. We investigate the dynamics of spike volleys with and without ephaptic coupling, and the resulting differences in axonal delays. There are several structural parameters that we keep fixed for simplicity, such as the fibre volume fraction (80% [36]), the fibre length (10cm), and the distribution of axon diameters. The spikes are generated at the proximal end of the fibre bundle, with spike times drawn from a uniform distribution. The width of this distribution determines the duration of a stimulus, and the number of spikes determines its intensity. We record the arrival time when a spike reaches the distal end of the axon, and the difference between the arrival time and the time the spike was initiated at the proximal end constitutes the axonal delay.

We first investigated how axonal delays are affected by ephaptic coupling, and focused on the mean of the delay distribution. In the presence of ephaptic coupling, we observe a decrease of the mean axonal delays as the stimulus intensity is increased (solid lines in Fig 8). In the absence of such coupling, the mean axonal delays remain constant (dashed lines in Fig 8). The stimulus duration is set to either 1ms, 2ms and 3ms, and the bundle diameters are varied between 2mm and 8mm. The mean axonal delays drop nonlinearly with increasing intensity in the presence of ephaptic coupling, but remain unchanged in its absence (Fig 8A–8C). At full intensity (100%) and with ephaptic coupling, the mean axonal delays drop from 35ms to 20ms as the diameter of the fibre bundle is increased from 2mm to 8mm if the stimulus duration is 1ms (Fig 8A). At 2ms stimulus duration, the mean axonal delays drop from 36ms (unchanged) to 28ms with increasing diameter of the fibre bundle (Fig 8B). In other words, the mean axonal delay decreased by up to 40% at full stimulus intensity.

thumbnail
Fig 8. Increasing the stimulus intensity, i.e. the number of spikes in a volley, decreases axonal transmission times and the latency of stimulus response.

A-C: Mean axonal delay with ephaptic coupling (solid) and without ephaptic coupling (dashed) for A: 1ms, B: 2ms, and C: 3ms stimulus duration. D-F: Standard deviation from the mean of axonal delay with ephaptic coupling (solid) and without ephaptic coupling (dashed) for D: 1ms, E: 2ms, and F: 3ms stimulus duration. Mean and standard deviation are computed from the distribution of delay times (cf. Fig 6D). G-I: Latency from stimulus onset to first maximum in neural mass model at G: 1ms, H: 2ms, and I: 3ms stimulus duration. Lines (shaded areas) indicate mean (1σ confidence interval) across 5 simulations. Colours indicate different bundle diameters.

https://doi.org/10.1371/journal.pcbi.1007858.g008

Next, we explored how the standard deviation of axonal delays (a measure for its dispersion) behaved in the presence of ephaptic coupling. We found that its qualitative behaviour is different from the mean of the axonal delays (Fig 8D–8F), with an initial decrease and a subsequent increase in the standard deviation.

Finally, we incorporated the axon bundle into the Jansen-Rit neural mass model [37]. The arrival of each spike at the distal end generates a current that is injected into the neural mass model. The response latency is determined by the time difference between stimulus onset and the maximum response of the neural mass model. This results in increased latencies as the stimulus duration is increased. However, in the presence of ephaptic coupling, we observe again a decrease in latencies as the stimulus intensity is increased, whereas in the absence of ephaptic coupling the decrease is only marginal (Fig 8G–8I). Regardless of stimulus duration, at full stimulus intensity ephaptic coupling reduces the response latency by up to 8ms, which corresponds to a reduction by approximately 15%.

Discussion

The key finding of our study is that spike volleys generate EPs with sufficiently large amplitudes to modulate axonal delays. Specifically, the mean delay of a spike volley decreases as the number of spikes in the spike volley is increased. Therefore, our results suggest that varying the amplitude of a neuronal signal can adjust its delay. Using a neural mass model, we have demonstrated that the decrease of axonal delays translates into the decrease of response latencies as the stimulus intensity is increased.

We have calibrated the spike propagation model using a biophysical model, by comparing the velocity change of a spike volley within a fibre bundle composed of identical axons. Here we observed the opposite effect: In the presence of ephaptic coupling, the spike volley slowed down. This is in line with previous numerical studies which investigated ephaptic coupling effects between a small number of identical axons. There, ephaptic coupling led to a synchronisation of the spikes within a volley, and a concurrent deceleration of the spikes. The acceleration of spike volleys that we observe in fibre bundles with distributed axon diameters can be attributed to the dispersive effect, which results from the axon diameter distribution and prevents synchronisation of a spike volley. As we show in Fig 5, an axially distributed spike volley causes primarily depolarisation within the fibre bundle, which then results in the acceleration of the majority of spikes within the volley. Therefore, our results do not contradict previous studies, but generalise previous modelling approaches. Nevertheless, the spike propagation model that we devised here required several assumptions that we are going to discuss in more detail.

We computed the EP using the line approximation (i.e. the axon is assumed to be infinitely thin), which has been demonstrated to be very accurate [30]. We further assumed that the axon bundle is large, circular, homogeneous, and densely populated with axons. The latter is justified by electron micrography studies which suggest that only a small fraction of an axon bundle is made up of extracellular space [35, 38]. Since axonal membranes and the myelin sheaths have a much larger resistivity compared to the extracellular medium, electric currents can only pass through the extracellular medium. We assumed that the medium is homogeneous and that the effective conductivity of the fibre bundle is the conductivity of the extracellular fluid multiplied by the relative size of the extracellular space. This calls for more detailed simulations of the spread of EPs with spatial heterogeneity taken into account. For mathematical convenience, we chose the fibre bundles to be large with circular cross sections. Realistic fibre bundles are indeed large, but often show a more sheet-like morphology [39]. An open question is whether this morphology influences the effect of EPs within our framework (a recent study used coupled axons with FitzHugh-Nagumo dynamics to demonstrate ephaptic coupling effects in sheet-like bundles [20]).

Furthermore, we ignored possible effects due to the axonal microstructure. We assumed the axonal membrane to be smooth (effectively a homogenised axon [40]), and that therefore nodes of Ranvier are not relevant as point sources. This is certainly the case at large distances from the axon, as can be seen in Fig 2F. However, at close proximity such effects would be relevant, as the EP at a node of Ranvier can reach several hundreds of microvolts. It is unknown whether nodes of neighbouring axons are sufficiently aligned to affect action potential generation in such a manner. As oligodendrocytes can myelinate multiple axons [41, 42], it is conceivable that neighbouring axons show some degree of alignments, in which case it would be possible to observe ephaptic coupling effects in much smaller fibre bundles, provided the spike volleys are sufficiently synchronised. Furthermore, we have demonstrated that computing the EP in a fibre bundle can be very well approximated by Eq (2). This expression depends only on the membrane potential V(z) and a convolution thereof. As can be seen from Fig 2, the spatial profile of the membrane potential is quite smooth along the axon, which is due to the fact that the length of a characteristic spike is by at least one order of magnitude larger than the length of a myelinated segment.

To demonstrate the effect of EPs on axonal delays we used a strongly simplified model for spike propagation. This model assumes that the spike velocity resulting from the axon morphology is known, and that this velocity is perturbed by the EP. It does not contain possible compensation effects arising from Hodgkin-Huxley dynamics (i.e. subthreshold currents that repolarise the axonal membrane), and further studies are required using the Hodgkin-Huxley framework to confirm our results. While we have used a simple scenario to calibrate the spike propagation model, to reproduce all of our findings within the Hodgkin-Huxley framework would be computationally extremely expensive, and is therefore beyond the scope of the present study.

We have incorporated the axon bundle model into the Jansen-Rit neural mass model to build a model system for primary sensory information processing, and to investigate the relationship between stimulus intensity and response latency. Psychological experiments across different sensory modalities yield the same qualitative relationship, whereby the latency decreases with increasing intensity [2127]. Such experiments typically measure the delay between stimulus presentation and the first maximum or minimum of the neural response measured electrographically. To replicate this experimental design, we measured the time difference between stimulus onset, i.e. the start of the spike volley, to the first maximum in the response of the Jansen-Rit model. Interestingly, only the presence of ephaptic coupling could explain the latency-intensity relationship. We are aware that the Jansen-Rit model is a fairly simple representation of a cortical microcircuit, and that other nonlinear processes not taken into account in its derivation may also reduce the response latency with increased stimulus intensity, such as oscillation-mediated information transmission [43]. Nevertheless, our modelling approach suggests that ephaptic coupling effects play a role in neural responses.

While there is such implicit evidence, further experimental studies are necessary to test our hypothesis. The experimental design would be highly invasive, since the EPs drop rapidly with distance outside fibre bundles. Animal experiments have already demonstrated the possibility to record EPs within axonal fibre bundles [32, 44]. An interesting test bed could also be a delay analysis within stimulation-response paradigms used in epileptic patients to determine the seizure focus [45, 46].

If such activity-dependent (or rate-dependent) delays occur in fibre bundles, then one may speculate as to their putative role in information processing. Since axonal delays are in general quite small (about 30ms in a 10cm long fibre bundle), the main effect should be on fast oscillations. It is indeed tempting to propose that such variable delays may have an effect on long-range gamma synchronisation, and that synchronisation patterns can be flexibly switched by changes in the amplitude of the transmitted spike volleys. We have found that ephaptic coupling can decrease delays by up to 10ms, which would be one half of the period of a gamma cycle at 50Hz. It has been demonstrated that delays are critical in shaping the functional architecture of the brain [4749], and ephaptic modulation of such delays could therefore play a role in flexibly synchronising distant brain areas.

Methods

In this section, we describe the mathematical framework underlying our study. We use a detailed microscopic description of the interaction between axonal fibres, and a leading-edge approximation which reduces the computational effort, but retains key properties of the interaction. We investigate a fibre bundle in which axons are coupled by the EPs generated by spikes. We first show how to compute the EPs generated by single spikes and spike volleys, and then present the framework of the spike propagation model.

Single fibre

In an open fibre bundle, the EP is determined by currents entering and leaving an axon. The radial currents around an axon can be inferred from the spatial profile of an action potential [32] given by the cable equation: (6) with I(z, t) being the radial (trans-membrane) currents, V(z, t) the membrane potential, a the axon radius and σi the conductivity of the intracellular medium. The axial dimension is represented by z, and time by t. The EP, denoted by ϕ for the single spike, can then be computed from the radial current via: (7) with σe being the extracellular conductivity, and d the radial coordinate measuring the distance from the axon. Inserting Eq (6) into Eq (7) yields (8) This integral is the convolution of the curvature of the action potential profile with the kernel G(zz′) = ((zz′)2 + d2)−1/2. In general, this integral has to be evaluated numerically. In order to obtain an analytical solution, we approximate the shape of an action potential by piecewise linear or piecewise quadratic functions.

Piecewise approximation of action potential profile.

In general the profile of an action potential has to be determined either numerically, or using spike-diffuse-spike formalisms. In the former case it is impossible to parameterise the profile, and in the latter the analytical expressions are still prohibitive to follow through with the calculations of the EP. Therefore, we present a formalism which approximates the profile of an action potential with either piecewise linear or piecewise quadratic functions. This method can be extended to arbitrary polynomial expressions, and is similar to curve-fitting with splines.

piecewise linear approximation

The simplest approximation of an action potential is given by two linear functions on two consecutive intervals, describing the rising and the falling phase of the action potential, respectively: (9) The first derivative of this approximation is piecewise constant with discontinuities at z = z0, z = z1 and z = z2. The second derivative is therefore (10) It is then straightforward to compute the EP: (11) (12) piecewise quadratic approximation

For the piecewise quadratic approximation, we divide the AP profile into three segments: (13)

Given z0, z1, z2, z3 and Vmax, there are four unknowns a1, a2, a3 and zmax. To ensure a smooth profile, we impose boundary conditions that assume V(z) is smoothly differentiable, i.e. , , , and . After some manipulation, we obtain: (14) (15) (16) (17)

The second derivative of the spatial profile is piecewise constant: (18)

The EP is then found to be (19) (20) (21)

Fibre bundle

Bundle with identical axons.

To compute an upper boundary of the EP produced by multiple action potentials, we assume that a perfectly synchronous spike-volley travels through a dense fibre bundle. All axons in this fibre bundle have the same diameter and are arranged in concentric rings (Fig 3A). At the centre of an empty grid position we compute the EP by summing ϕ at distance (2n + 1)a of 6n axons, with n ranging from 1 to N, with N large: (22) Although ϕ is approximately 20μV at the surface of an isolated axon, in a fibre bundle the combined effect can lead to EPs of many mV. Interestingly, we find that for large enough fibre bundle diameters the profile of the cumulative EP is almost proportional to the profile of the generating action potentials. We give a mathematical explanation for this next.

Analytical solution.

The cumulative EP at the core of an axon bundle is computed with the following integral, (23) with P being the axon bundle diameter, and Ω(a) being the cross-sectional area occupied by an axon with radius a. We set Ω(a) = πa2/(ρg2), where g is the g-ratio and ρ is the fibre volume fraction.

Inserting Eq (1) into Eq (23), and solving the integral over ρ, results in (24)

Integration by parts then yields (25)

Next, we use the approximation (26) which leads to (27)

Using integration by parts, this integral ultimately yields (28)

We may regard this result as the far-field approximation of EPs in axonal fibre bundles.

We note here also that in the limit P → 0, exp(|zz′|/P)/P → 2δ(zz′), which yields (29)

Eq (24) suggests that the far-field approximation of the cumulative EP is independent of the axon morphology. At this point, however, we have not taken into account that axons of different diameters transmit action potentials at different velocities, and that therefore the spatial profile widens with increasing action potential velocity.

Spike propagation model

The two most common ways to model axonal signal transmission are either Hodgkin-Huxley type dynamics embedded in a core-conductor model, or simpler spike-diffuse-spike approaches. Both ways allow one to determine the spike velocity as a function of electrophysiological and structural parameters. Here, we employ a much simpler model that describes the position zi of an action potential (more precisely, its leading edge or rising phase) travelling along the ith axon by one simple equation: (30) where vi(z, t) is the velocity of the action potential as function of the axial direction z and time t. If the axon is homogeneous and does not experience spatial or temporal perturbations, then the velocity can be expressed by vi(z, t) = vi,0, which is the intrinsic velocity of the axon, determined by its morphological and electrophysiological properties. We assume here that this velocity is known for each axon. In the absence of perturbations, one can therefore express the axonal delays by τi = L/vi,0, with L being the length of the fibre bundle. We set here vi,0 = αdi, with di being the diameter of the ith axon, and α = 5ms−1/μm.

Changes in the EP lead to perturbations of the membrane potential of an axon. A negative (positive) EP effectively depolarises (hyperpolarises) the axonal membrane, and therefore increases (decreases) the propagation velocity. A convenient formalism to incorporate such changes is the spike-diffuse-spike framework, in which the spiking threshold is a parameter to explicitly describe the onset of an action potential [19]. Such thresholds can also be determined within the Hodgkin-Huxley formalism, albeit these thresholds vary with the depolarisation rate [50]. The EP can therefore be regarded as a perturbation of such a threshold. Within the spike-diffuse-spike framework, the following relationship between the spike threshold Vthr and the delay of spike generation Δt between two consecutive nodes of Ranvier can be derived: (31) see Fig 6F for a visual representation. The function ft) depends on structural and electrophysiological parameters of the axon. The EP can be incorporated into the spiking threshold, Vthr(z, t) = Vthr,0 + EP(z, t), with Vthr,0 being the uniform spiking threshold of the unperturbed axon. Via Eq (31) one can relate Vthr,0 to Δt0 of the unperturbed axon, and to its intrinsic velocity via vi,0 = lt0, where l is the distance between two consecutive nodes of Ranvier. To further simplify our scheme we linearise Eq (31) around Vthr,0 and Δt0: (32) which can be reformulated into (33) Here we made use of v = Δzt(z, t), which results in Δt(z, t)/Δt0 = v0/v(z, t). The parameter c denotes the relative steepness of ft) around Δt0. For simplicity, we lump c and Vthr into one parameter γ = 1/(cVthr), which results in Eq (3). This parameter can be used as a tuning parameter for the inverse strength of the ephaptic coupling, and in the numerical simulations we set γ = 0 to represent the absence of any ephaptic interaction.

The computation of the EP generated by a spike requires knowledge about its axial profile, which can be inferred from the temporal evolution of a spike at any given point along the axon. In order to obtain a computationally efficient framework, we focus on the piecewise linear approximation. The temporal profile of a spike in linear approximation is given by (34) where t0 is the time of onset of the spike, t1 is the time it reaches its maximum depolarisation, and t2 is the time when the membrane is fully repolarised. The variable v(z, t) represents the velocity of the leading edge of the spike, which can vary rapidly. We introduce the effective velocity to estimate the velocity of the entire spike, or of its centre of mass. This effective velocity is, in essence, a time-averaged quantity. We approximate the effective velocity of a spike by (35) with τ = 1ms. The temporal profile can then be translated into the spatial profile by setting z = veff(z) × t, and zm = veff(z) × tm. The parameters for the spike propagation model are listed in Table 1.

thumbnail
Table 1. List of parameters used for the spike propagation model.

https://doi.org/10.1371/journal.pcbi.1007858.t001

Biophysical model

An axon can be regarded as a core-conductor with intracellular potential ϕi(z, t), and extracellular potential ϕe(z, t). The difference between in intracellular and extracellular potential results in the membrane potential V(z, t) = ϕi(z, t) − ϕe(z, t). The extracellular potential is subject to the model setup. Here we consider an axon bundle composed of identical axons, a fraction of which carries spikes synchronously. Because the extracellular potential is considered homogeneous within the fibre bundle (Eq 2), the dynamics of all spikes is identical, and it suffices to solve one cable equation to describe the propagation of all spikes: (36) Since we consider myelinated axons, the radial capacitance Cm and the radial resistance Rm vary with the axial coordinate z, depending on whether a segment is myelinated (internode) or unmyelinated (node of Ranvier): (37) The model axon is myelinated periodically with nodes of length l, capacitance Cn, and resistivity Rn, and myelinated internodes with length L, capacitance Cmy and resistivity Rmy. The intracellular resistance Rax remains constant along the axon.

The term I(V) represents voltage-gated currents modelled with the Hodgkin-Huxley model. These currents only occur in nodal segments: (38) The gating variables m, h, p, s obey the following dynamics: (39) with n = m, h, p, s. The variables αn and βn are defined as follows:

The parameters for the biophysical model are given in Table 2.

thumbnail
Table 2. List of parameters used for the biophysical model.

https://doi.org/10.1371/journal.pcbi.1007858.t002

Jansen-Rit microcircuit

The spike volleys represent cortical, subcortical or sensory information being transitted by axon bundles. To describe the response of neuronal circuits, e.g. cortical microcircuits, we use the Jansen-Rit model [37, 51, 52] and record the maximum response in the membrane potential of its pyramidal cell population. The Jansen-Rit model is composed of six differential equations: (40) Here, y0 is the postsynaptic potential (PSP) generated by the output of the pyramidal cells at the two interneuron types, and y1 and y2 are the excitatory and inhibitory PSPs generated at the pyramidal cells by external stimuli p(t) and the firing activity of the interneurons. Furthermore, y3, y4 and y5 are auxiliary variables in the synaptic conversion of firing rates into PSPs, with a and b being the inverse time constants of excitatory and inhibitory synapses, and A and B the respective maximum amplitudes of the synaptic response. The PSPs are converted into firing rates by the sigmoidal function (41) with e0 being the maximum firing rate, v0 the membrane potential at half of the maximum firing rate, and r sets the steepness of the sigmoid. The interaction between the different neuron types is scaled by the connectivity constants C1 to C4. The parameters are chosen as in [37] (see Table 3).

thumbnail
Table 3. List of parameters used for the Jansen-Rit model.

https://doi.org/10.1371/journal.pcbi.1007858.t003

The external firing rate p(t) is generated by the incoming spikes, (42) with tn being the arrival time of the nth spike, N the total number of spikes, and P = 0.1s−1 sets the synaptic coupling strength for the spike train.

The membrane potential y of the pyramidal cell population is determined by the difference between excitatory and inhibitory PSPs, i.e. y = y1y2. The response latency is then calculated as the position of the absolute maximum of y.

References

  1. 1. Buzsáki G, Anastassiou CA, Koch C. The origin of extracellular fields and currents—EEG, ECoG, LFP and spikes. Nat Rev Neurosci. 2012;13:407–420. pmid:22595786
  2. 2. Arvanitaki A. Effects evoked in an axon by the activity of a contiguous one. J Neurophysiol. 1942;5:89–108.
  3. 3. Anastassiou CA, Perin R, Markram H, Koch C. Ephaptic coupling of cortical neurons. Nat Neurosci. 2011;14:217–223. pmid:21240273
  4. 4. Anastassiou CA, Koch C. Ephaptic coupling to endogenous electric field activity: why bother? Curr Opin Neurobiol. 2015;31:95–103. pmid:25265066
  5. 5. Destexhe A, Contreras D, Steriade M. Spatiotemporal analysis of local field potentials and unit discharges in cat cerebral cortex during natural wake and sleep states. J Neurosci. 1999;19:4595–4608. pmid:10341257
  6. 6. Lindén H, Tetzlaff T, Potjans TC, Pettersen KH, Grün S, Diesmann M, et al. Modeling the spatial reach of the LFP. Neuron. 2011;72:859–872. pmid:22153380
  7. 7. Łȩski S, Lindén H, Tetzlaff T, Pettersen KH, Einevoll G. Frequency dependence of signal power and spatial reach of the local field potential. PLoS Comp Biol. 2013;9:e1003137. pmid:23874180
  8. 8. Einevoll G, Kayser C, Logothetis NK, Panzeri S. Modelling and analysis of local field potentials for studying the function of cortical circuits. Nat Rev Neurosci. 2013;14:770–785. pmid:24135696
  9. 9. Herreras O. Local Field Potentials: Myths and Misunderstandings. Front Neural Circuits. 2016;10:101. pmid:28018180
  10. 10. Katz B, Schmitt OH. Electric interaction between two adjacent nerve fibres. J Physiol. 1940;97:471–488. pmid:16995178
  11. 11. Rosenblueth A. The stimulation of myelinated axons by nerve impulses in adjacent myelinated axons. Amer J Physiol. 1941;132:119–128.
  12. 12. Marrazzi AS, Lorente de Nó R. Interaction of neighbouring fibres in myelinated nerve. J Neurophysiol. 1944;7:83.
  13. 13. Eilbeck JC, Luzader SD, Scott AC. Pulse evolution on coupled nerve fibres. Bull Math Biol. 1981;43:389–400. pmid:7317672
  14. 14. Barr RG, Plonsey R. Electrophysiological interaction through the interstitial space between adjacent unmyelinated fibers. Biophys J. 1992;61:1164–1175. pmid:1600078
  15. 15. Binczak S, Eilbeck JC, Scott AC. Ephaptic coupling of myelinated nerve fibers. Physica D. 2001;148:159–174.
  16. 16. Reutskiy S, Rossoni E, Tirozzi R. Conduction in bundles of demyelinated nerve fibers: computer simulation. Biol Cybern. 2003;89:439–448. pmid:14673655
  17. 17. Maına I, Tabi CB, Ekobena Fouda HP, Mohamadou A, Kofané TC. Discrete impulses in ephaptically coupled nerve fibers. Chaos. 2015;25:043118. pmid:25933666
  18. 18. Goldwyn JH, Rinzel J. Neuronal coupling by endogenous electric fields: cable theory and applications to coincidence detector neurons in the auditory brain stem. J Neurophysiol. 2016;115:2033–2051. pmid:26823512
  19. 19. Schmidt H, Knösche TR. Action potential propagation and synchronisation in myelinated axons. PLoS Comp Biol. 2019;15:e1007004. pmid:31622338
  20. 20. Sheheitli H, Jirsa VK. A mathematical model of ephaptic interactions in neuronal fiber pathways: Could there be more than transmission along the tracts? Network Neurosci. 2020;4:595–610. pmid:32885117
  21. 21. Bell AH, Meredith MA, Van Opstal AJ, Munoz DP. Stimulus intensity modifies saccadic reaction time and visual response latency in the superior colliculus. Exp Brain Res. 2006;174:53–59. pmid:16528494
  22. 22. Adler G, Adler J. Auditory Stimulus Processing at Different Stimulus Intensities as Reflected by Auditory Evoked Potentials. Biol Psychiatry. 1991;29:347–356. pmid:2036478
  23. 23. Jaskownki P, Rybarczyk K, Jaroszyk F. The relationship between latency of auditory evoked potentials, simple reaction time, and stimulus intensity. Psychol Res. 1994;56:59–65.
  24. 24. Ulrich R, Rinkenauer G, Miller J. Effects of Stimulus Duration and Intensity on Simple Reaction Time and Response Force. J Exp Psychol. 1998;24:915–928. pmid:9627425
  25. 25. Bergamin O, Kardon RH. Latency of the Pupil Light Reflex: Sample Rate, Stimulus Intensity, and Variation in Normal Subjects. Invest Ophthalmol Vis Sci. 2003;44:1546–1554. pmid:12657591
  26. 26. Chapman CR, Oka S, Bradshaw DH, Jacobson RC, Donaldson GW. Phasic pupil dilation response to noxious stimulation in normal volunteers: Relationship to brain evoked potentials and pain report. Psychophysiology. 1999;36:44–52. pmid:10098379
  27. 27. Mayhew SD, Iannetti GD, Woolrich MW, Wise RG. Automated single-trial measurement of amplitude and latency of laser-evoked potentials (LEPs) using multiple linear regression. Clin Neurophys. 2006;117:1331–1344. pmid:16644270
  28. 28. Diesmann M, Gewaltig MO, Rotter S, Aertsen A. State space analysis of synchronous spiking in cortical neural networks. Neurocomputing. 2001;38–40:565–571.
  29. 29. Moldakarimov S, Bazhenov M, Sejnowski TJ. Feedback stabilizes propagation of synchronous spiking in cortical neural networks. Proc Nat Ac Sci. 2007;112:2545–2550.
  30. 30. Holt GR, Koch C. Electrical Interactions via the Extracellular Potential Near Cell Bodies. J Comp Neurosc. 1999;6:169–184. pmid:10333161
  31. 31. Gold C, Henze DA, Koch C, Buzsáki G. On the origin of the extracellular action potential waveform: A modelling study. J Neurophysiol. 2006;95:3113–3128. pmid:16467426
  32. 32. McColgan T, Liu J, Kuokkanen PT, Carr CE, Wagner H, Kempter R. Dipolar extracellular potentials generated by axonal projections. eLife. 2017;6:e26106. pmid:28871959
  33. 33. Arancibia-Cárcamo IL, Ford MC, Cossell L, Ishida K, Tohyama K, Atwell D. Node of Ranvier length as a potential regulator of myelinated axon conduction speed. eLife. 2017;6:e23329. pmid:28130923
  34. 34. Phillips KA, Stimpson CD, Smaers JB, Raghanti MA, Jacobs B, Popratiloff A, et al. The corpus callosum in primates: processing speed of axons and the evolution of hemispheric asymmetry. Proc R Soc B. 2015;282:20151535. pmid:26511047
  35. 35. Liewald D, Miller R, Logothetis N, Wagner HJ, Schüz A. Distribution of axon diameters in cortical white matter: an electron-microscopic study on three human brains and a macaque. Biol Cybern. 2014;108:541–557. pmid:25142940
  36. 36. Xu J, Li H, Harkins KD, Jiang X, Xie J, Kang H, et al. Mapping mean axon diameter and axonal volume fraction by MRI using temporal diffusion spectroscopy. Neuroimage. 2014;103:10–19. pmid:25225002
  37. 37. Jansen BH, Rit VG. Electroencephalogram and visual evoked potential generation in a mathematical model of coupled cortical columns. Biol Cybern. 1995;73:357–366. pmid:7578475
  38. 38. Aboitiz F, Scheibel AB, Fisher RS, Zaidel E. Fiber composition of the human corpus callosum. Brain Res. 1992;598:143–153.
  39. 39. Wedeen VJ, Rosene DL, Wang R, Dai G, Mortazavi F, Hagmann P, et al. The Geometric Structure of the Brain Fiber Pathways. Science. 2012;335:1628–1634. pmid:22461612
  40. 40. Basser PJ. Cable equation for a myelinated axon derived from its microstructure. Med & Biol Eng Comput. 1993;31:S87–S92. pmid:8231331
  41. 41. Simons M, Nave KA. Oligodendrocytes: Myelination and Axonal Support. Cold Spring Harb Perspect Biol. 2016;8:a020479.
  42. 42. Harty BL, Coelho F, Pease-Raissi SE, Mogha A, Ackerman SD, Herbert AL, et al. Myelinating Schwann cells ensheath multiple axons in the absence of E3 ligase component Fbxw7. Nat Commun. 2019;10:2976. pmid:31278268
  43. 43. Hahn G, Ponce-Alvarez A, Deco G, Aertsen A, Kumar A. Portraits of communication in neuronal networks. Nat Rev Neurosci. 2019;20:117–127. pmid:30552403
  44. 44. Kuokkanen PT, Ashida G, Kraemer A, McColgan T, Funabiki T, Wagner H, et al. Contribution of action potentials to the extracellular field potential in the nucleus laminaris of barn owl. J Neurophysiol. 2018;119:1422–1436. pmid:29357463
  45. 45. Matsumoto R, Nair DR, LaPresto E, Bingaman W, Shibasaki H, Lüders HO. Functional connectivity in human cortical motor system: a cortico-cortical evoked potential study. Brain. 2007;130:181–197. pmid:17046857
  46. 46. Trebaul L, Deman P, Tuyisenge V, Jedynak M, Hugues E, Rudrauf D, et al. Probabilistic functional tractography of the human cortex revisited. Neuroimage. 2018;181:414–429. pmid:30025851
  47. 47. Petkoski S, Spiegler A, Proix T, Aram P, Temprado JJ, Jirsa VK. Heterogeneity of time delays determines synchronization of coupled oscillators. Phys Rev E. 2016;94:012209. pmid:27575125
  48. 48. Petkoski S, Jirsa VK. Transmission time delays organize the brain network synchronization. Phil Trans R Soc A. 2019;377:20180132. pmid:31329065
  49. 49. Finger H, Gast R, Gerloff C, Engel AK, König P. Probing neural networks for dynamic switches of communication pathways. PLoS Comp Biol. 2019;15:e1007551. pmid:31841504
  50. 50. Platkiewicz J, Brette R. A Threshold Equation for Action Potential Initiation. PLoS Comp Biol. 2010;6:e1000850. pmid:20628619
  51. 51. Grimbert F, Faugeras O. Bifurcation analysis of Jansen’s neural mass model. Neural Comput. 2006;18:3052–3068. pmid:17052158
  52. 52. Spiegler A, Knösche TR, Schwab K, Haueisen J, Atay FM. Modelling Brain Resonance Phenomena Using a Neural Mass Model. PLoS Comp Biol. 2011;7:e1002298.