Fields, Geometry and Their Impact on Axon Interaction

Aman Chawla, Salvatore Domenic Morgera, Arthur D. Snider

Alliance University, India and University of South Florida, Tampa, USA.

DOI: 10.4236/jamp.2021.94053                

Abstract

In this paper, we present a simulation program that allows for the concurrent propagation of action potentials in axons coupled via currents, as well as, for the first time, the computation of the resultant nodal electric field generated as the action potentials traverse the tract of axons. With these fields in hand, we inject currents into nodes of axons that depend on these fields and study the coupling between axons in the presence of the fields and currents present jointly in varying strengths. We find close-to-synchronized propagation in three dimensions. Further, we derive for the first time a mathematical equation for tortuous tracts (as opposed to linear) with such field-mediated coupling. The geometrical formulation enables us to consider spacetime perturbative effects, which have also not been considered in the literature so far. We investigate the case when gravitational radiation is present, in order to determine its impact on tract information processing. We find that action potential relative-timing in a tract is affected by the strength and frequency of gravitational waves and the waning of this influence with weakening strength. This latter study blurs the division between what lies inside and outside man. As an additional novelty, we investigate the influence of geometry on the information transmission capacity of the ephaptically-coupled tract, when viewed as a discrete memoryless channel, and find a rising trend in capacity with increasing axonal inclinations, which may occur in traumatic CNS injury.

Keywords

Field-Mediated Axon Interaction, Current-Mediated Axon Interaction, Metric Perturbations and Axon Interaction, Capacity of Interacting Axons

1. Introduction

What would your consciousness (or thoughts, more simply) be like when approaching a black hole? Or perhaps when your brain is subject to strong electromagnetic fields? Can we tackle such questions scientifically? If not the subjective perception itself, then perhaps we may enquire about just the brain state in general in such extreme conditions. The premise of this paper is that if we can study nerve tracts from the perspective of the four forces (electromagnetic, strong, weak and gravitational), then we can perhaps lay our finger on and pin-point exactly what makes these tracts enablers of consciousness.

Though the mind may be more than just a brain state [1], when making a connection between brain states and consciousness we rely on the observation that brain states are nevertheless a window to the mind [2]. This observation was made in the context of EEG signals which are electrical in nature, as being representative of the underlying brain state. But consider also the fact that EEG signals are more than just physical signatures—they carry bits of information via communication signals. Neither view, the EEG-physical view nor the information-communication view precludes the other. Viewed one way, via an electric voltage meter, the brain is generating electric voltages. At this level, there is no hint of “bits”. However, when we superpose the picture that certain voltage waveforms “mean” particular bit values, then we see bits in all neural signals. In Minskian terminology, we have gone from an `agent’ perspective (the physical voltages) to an “agency” perspective where the bits can now participate in information processing which subsequently directs various activities of the organism.

The Minskian agent-agency theory is not the only view that can be taken in this regard. For example, Landauer [3] uses entirely different arguments to convince the reader that bits are physical, establishing a sort of identity or closer relation between the representation (bits) and the representative (voltages of the physical substrate). Regardless of whether one prefers the Minskian view or the Landauer perspective, bits are important and thus in this paper we also include information as one of the important modalities to be investigated. Further, this modality is also important since there is a line of philosophical thinking which equates “me” with my brain’s state of information processing [4].

The setting in which our work is carried out is as follows. We know that the brain is composed primarily of neurons and glial cells. A neuron is often seen as a tri-partite structure, with dendritic arbor, cell body and axon. Suppose two axons are considered, placed in extracellular space. Action potentials may be computed and would hop from node (of Ranvier) to node on each of the axons, using the mechanism of saltatory conduction [5], with the action potentials interacting with each other [6]. Concomitant with these action potentials, electric fields would be generated at each node of Ranvier. All these fields can be computed and superposed with one another in order to determine the net electrodynamics of the tract. In this work, we take our reader through a few accurate simulations of such action potential propagation in axon tracts in the presence of the concomitant internally generated electric fields (related to Local Field Potentials (LFPs)), which has many implications for improving our understanding of neuroscience. It also has several applications related to the correct assessment of diseased conditions of the nervous system, such as multiple sclerosis (MS) and autism spectrum disorder (ASD).

Along the way, we also look at impinging gravitational radiation and issues of information processing (channel capacity, in particular) in the context of axon tracts, all from the computational angle—bolstered by numerical simulations. While there is extant literature on electric and magnetic fields in volume conductors [2], the literature on gravitational waves interacting with biophysical systems is rather sparse—or non-existent. Here we will only briefly comment on the relationship of our work to [7]. That work is concerned with the representation of spacetime in the brain, whereas we are concerned first of all with the actual impact of spacetime perturbations on information processing by the three-dimensional axon tract. However, if we model our three-dimensional information processing along the lines suggested in that paper, by connecting internal firing activity patterns further to external spacetime reality, then we can build a bridge between that work and our present work. This would be interesting since we would be looking at the impact of spacetime perturbations on the very perception of that spacetime itself.

This paper is structured as follows. In the next section, Section 2, we look at electric field-mediated interaction between axons. In Section 3, we review current-mediated interaction. In Section 4, we introduce the new formalism for tract tortuosity. In Section 5, we simulate the field-mediated interaction due to sodium effects alone and in Section 6, due to both sodium and potassium ions. In Section 7, we turn to gravitational radiation’s impact and in Section 8, we look at issues of tract capacity. Finally, we conclude in Section 9, providing directions for future investigation.

Notational Preliminaries

We define our notation in Table 1 and Table 2.

2. Axon-Axon Interaction via Electric Fields

In this introductory section, we will study the interaction between axons as mediated by the electric fields generated by the very axons themselves. This is therefore a sort of feedback coupling within the axon tract, and, as our results indicate, if the feedback is included positively (by means of a gain which will be discussed later), it leads to a dramatic speed up of the action potential propagation. For background on much of this section, including a few of the results, we refer the reader to [8].

2.1. Nature of the Current Source at a Single Ion Channel

We are interested in the current conducted through an ion channel. Though it will depend on the voltage across the ion channel, it is useful to look at the velocity of ions in the channel to get an idea as to whether these ions approach high speeds or not. From [9], we know that as per classical laws (i.e. not invoking quantum mechanics), there are roughly 3.1 × 10⁷ ions per second flowing

Table 1. List of principal symbols and their meanings I.

Table 2. List of principal symbols and their meanings II.

through a typical sodium ion channel, which corresponds to a sodium current of about 5 pico-Amperes. This is close to the value of 2 pico-Amperes as is commonly measured using patch clamp techniques. This also means that a single ion takes about 32.25 nano-seconds to cross an ion channel pore. The length of such a pore is about 5 Angstroms or 0.5 nano-meters. This leads to a mean velocity through the pore of $v_{ion}=\frac{0.5\times {10}^{-9}}{32.25\times {10}^{-9}}\frac{\text{m}}{\text{s}}=15.5\times {10}^{-3}\frac{\text{m}}{\text{s}}$, placing the velocity firmly in the non-relativistic regime.

2.2. The Electric near Field Generated by a Sodium Ion Channel Ring

In this section, our goal will be to compute the electric field that results from an open sodium ion channel, which permits a peak current of 2 pico-Amperes through itself. Since axon fibers are tightly packed in a tract, with very small mutual separation, what we are interested in may be considered to be the near-field generated by a node of Ranvier. Note that the term near-field is usually reserved for radiating antennas whose largest dimension and radiation wavelength can be used to compute the boundary of the near field zone, via the Fraunhofer distance [10].

For the case at hand the near field boundary is about 2 to 3 microns from the surface of the node of Ranvier (see page e4 in [11] for the refractive index of brain tissue which can be used to easily compute the Fraunhofer distance). If axons are tightly packed, this would imply that they are all within each other’s near field. We develop our results for the sodium ion¹. The next equation, an ansatz based on volume conduction theory (as presented in Chapter 2.2 of [12] ), models the volume source current density produced by an ion channel ring (Figure 1), in cylindrical coordinates.

$I_v\left(\stackrel{\to }{r^{\prime }},t\right)=I_0\left(1-\frac{{\rho }^{\prime }}{p}\right)\delta \left(z^{\prime }\right)U\left(t-r_1\right)U\left(l_1-t\right){\hat{\rho }}^{\prime }$ (1)

Next, to find the electric near field associated with this current density, we plug it into [13]

${\varphi }_{LFP}\left(\stackrel{\to }{r},t\right)=\frac{1}{4\pi \epsilon }{\displaystyle \int }\frac{I_v\left(\stackrel{\to }{r^{\prime }},t\right)}{\left|\stackrel{\to }{r}-\stackrel{\to }{r^{\prime }}\right|}{\text{d}}^3{r}^{\prime }$ (2)

which gives the LFP as a volume integral involving the current density of the source and the distance between the source and the field positions. Next, using the quasi-static approximation, $\stackrel{\to }{E}=-\nabla \varphi$ gives us (see [14] [15] ) the electric near field as the negative gradient of the potential computed in Equation (2). To justify the use of the previous two equations, we quote [16]:

A simple solution is to realize that in the near-field world, where the structures are much smaller than the wavelength, the electric field is mostly unrelated to the velocity of light; instead, the dominant component of the electric field is longitudinal, i.e., the quasistatic gradient of the Coulomb integral over the instantaneous charge distribution.

Using vector analysis [17] we have,

$\stackrel{\to }{E_{\text{channel}}}\left(\stackrel{\to }{r}\right)=\hat{k}\frac{z}{4\pi \epsilon }{\displaystyle \int }\frac{I\left(\stackrel{\to }{r^{\prime }}\right)}{{\left|\stackrel{\to }{r}-\stackrel{\to }{r^{\prime }}\right|}^3}{\text{d}}^3\stackrel{\to }{r^{\prime }}+\cdots$ (3)

${\hat{e}}_{\rho }\frac{1}{4\pi \epsilon }{\displaystyle \int }\frac{I\left(\stackrel{\to }{r^{\prime }}\right)\left(\rho -{\rho }^{\prime }\cos \left(\varphi -{\varphi }^{\prime }\right)\right){\text{d}}^3\stackrel{\to }{r^{\prime }}}{{\left|\stackrel{\to }{r}-\stackrel{\to }{r^{\prime }}\right|}^3}$ (4)

By shifting the ansatz, Equation (1), outward from the origin by an amount ${{\rho }^{\prime }}_1$, we can write:

Figure 1. Ion channel ring at a node of Ranvier of an axon. The indentations are the ion channels, placed around the ring at various points.

$I\left(\stackrel{\to }{r^{\prime }}\right)=(\begin{array}{c}\begin{array}{l}00<{\rho }^{\prime }<{{\rho }^{\prime }}_1\hfill \\ \delta \left(z^{\prime }\right)\frac{I_0}{{{\rho }^{\prime }}_2-{{\rho }^{\prime }}_1}\left({{\rho }^{\prime }}_2-{\rho }^{\prime }\right){{\rho }^{\prime }}_1<{\rho }^{\prime }<{{\rho }^{\prime }}_2\hfill \\ 0{\rho }^{\prime }>{{\rho }^{\prime }}_2\hfill \end{array}\end{array}$ (5)

Next, we note the following standard result on the integral of a Dirac delta function:

${\displaystyle \underset{-\infty }{\overset{\infty }{\int }}\delta \left(z^{\prime }\right)\text{d}{z}^{\prime }}=\mathrm{1,}$ (6)

a standard result from basic vector calculus:

${\text{d}}^3\stackrel{\to }{r^{\prime }}={\rho }^{\prime }\text{d}{\rho }^{\prime }\text{d}{\varphi }^{\prime }\text{d}z\mathrm{,}$ (7)

and an application of the cosine formula from trigonometry:

$\left|r-r^{\prime }\right|=\sqrt{{\rho }^2+{{\rho }^{\prime }}^2-2\rho {\rho }^{\prime }\cos \left(\varphi -{\varphi }^{\prime }\right)+z^2}\mathrm{,}$ (8)

which, when used in the first term in Equations (3) and (4) yield

$\frac{I_0}{{{\rho }^{\prime }}_2-{{\rho }^{\prime }}_1}{\displaystyle \underset{{\varphi }^{\prime }=0}{\overset{\Delta {\varphi }^{\prime }}{\int }}{\displaystyle \underset{{\rho }^{\prime }={{\rho }^{\prime }}_1}{\overset{{{\rho }^{\prime }}_2}{\int }}\frac{\left({{\rho }^{\prime }}_2-{\rho }^{\prime }\right){\rho }^{\prime }\text{d}{\rho }^{\prime }\text{d}{\varphi }^{\prime }}{{\left[{{\rho }^{\prime }}^2+{\rho }^2-2\rho {\rho }^{\prime }\cos \left(\varphi -{\varphi }^{\prime }\right)+z^2\right]}^{\frac{3}{2}}}}}.$ (9)

Using a similar derivation, the second term in Equations (3) and (4) becomes

$\frac{I_0}{{{\rho }^{\prime }}_2-{{\rho }^{\prime }}_1}{\displaystyle \underset{{\varphi }^{\prime }=0}{\overset{\Delta {\varphi }^{\prime }}{\int }}{\displaystyle \underset{{\rho }^{\prime }={{\rho }^{\prime }}_1}{\overset{{{\rho }^{\prime }}_2}{\int }}\frac{\left({{\rho }^{\prime }}_2-{\rho }^{\prime }\right)\left(\rho -{\rho }^{\prime }\cos \left(\varphi -{\varphi }^{\prime }\right)\right){\rho }^{\prime }\text{d}{\rho }^{\prime }\text{d}{\varphi }^{\prime }}{{\left[{{\rho }^{\prime }}^2+{\rho }^2-2\rho {\rho }^{\prime }\cos \left(\varphi -{\varphi }^{\prime }\right)+z^2\right]}^{\frac{3}{2}}}}}.$ (10)

These two integrals are evaluated at user-input values of $\rho$, $\varphi$ and z. Using these evaluations, we compute the electric field, which then leads to coupling as will be discussed in Section 4. The more rudimentary means of axon-axon coupling is via extracellular currents, and this is reviewed next.

3. Axon-Axon Interaction via Currents: A Quick Pass

We begin by noting down the differential equation that we need to simulate, examining it, and giving a few remarks on its features. This equation is displayed here [6]:

$c_{m,k}\left(x\right)\frac{\partial V_k\left(x,t\right)}{\partial t}=\frac{1}{r_f}\frac{{\partial }^2{V}_k\left(x,t\right)}{\partial x^2}-\frac{{\alpha }_N}{r_f}{\displaystyle \underset{n=1}{\overset{N}{\sum }}}\frac{{\partial }^2{V}_n\left(x,t\right)}{\partial x^2}-\cdots$ (11)

$g_{m\mathrm{,}k}\left(x\right)V_k\left(x\mathrm{,}t\right)$ (12)

where $k=1,\cdots ,N$. The equation contains partial derivatives of the voltage in space and time. It also contains a number of constants. This system can be specifically classified as a system of coupled second-order linear partial differential equations (PDEs) in two independent variables and of the parabolic type, though simple parabolic equations do not contain terms like the last term on the right hand side of Equation (11) and Equation (12); that term is akin to a forcing term for the equation system [18].

In order to derive this equation, the authors of [6] first set up a current-balance between the intracellular and extracellular currents. Combining this with transmembrane potential differences, gives them a term for the ephaptic current. They include this ephaptic current term in the Frankenhaeuser-Huxley equation [19] for the nodal regions as well as the internodal regions. These two regions have two separate equations, each of which is the governing law for a specific set of longitudinal positions along each axon. The nodal equation contains the machinery for action potential generation via the ionic current term. External stimulation can also be accepted at a node. We refer the reader to [6] for a full understanding of the derivation of this equation.

In [20], we extend Equations (11) and (12) to the three-dimensional setting. We consider arbitrary inclinations for axons in a 3D coordinate system, but the axons are still linear themselves. We also consider arbitrary inter-axonal distances encoded in a specific matrix, which we label as W. We indicate how the new equation² is a generalization of the Frankenhaeuser-Huxley equation as well as Equations (11) and (12). In the next section, we embellish this equation to allow for tortuous axons and field-based interactions.

4. Generalized Ephaptic Equation with Tract Tortuosity and Electric Field-Based Coupling

>

One of our original goals with the three-dimensional axon tract formulation³ was to ultimately look at tortuous axons, axons wending their way along splines (Figure 2). Tortuosity of axons is a very important feature to be modeled as it substantially increases the length and surface area of axons. It has a purpose, not yet understood. In the generalized ephaptic equation (GEE) published in [20], if we make the angular inclinations and W-matrix z-dependent then that would take care of such tortuosity as well. The reasoning which leads us to this insight is that regardless of how tortuous a tract’s axons may be, if we cut the tract up into small segments of infinitesimal length dz the axons would be linear within those segments and have a fairly constant inclination as well. Then, stacking together these segments of length dz back to back, we would get a reconstruction

Figure 2. A spline is a smooth path, not necessarily linear. Axons quite often follow splines, and so modeling them as such is essential.

of the entire tortuous tract. Segment to segment, the inclination would gradually vary with z and so would the W-matrix entries. This leads us to write down the following equation for a tract of tortuous, coupled, axons.

$c_{m,k}\left(x\right)\frac{\partial V_k}{\partial t}=G_k^{\text{axoplasm}}{\cos }^2{\theta }_k\left(z\right)\left[\frac{{\partial }^2{V}_k}{\partial z^2}-{\displaystyle \underset{p=1}{\overset{N}{\sum }}}{W}_{ip}\left(z\right)K_p\left(z\right)\frac{{\partial }^2{V}_p}{\partial z^2}\right]-\cdots$ (13)

$G_{m\mathrm{,}k}{V}_k+I_k^{\text{injected}}$ (14)

wherein ${\theta }_k\left(z\right)$ is the angular inclination of axon k at the longitudinal position z along the tract axis, and $W_{ip}\left(z\right)$ is the interaxonal distance matrix at the same position. $K_p\left(z\right)$ depends on the angular inclination and is therefore also written as a function of z.

Next we modify the tortuous-tract equation to allow variable proportions of field-based and current-based coupling.

$c_{m,k}\left(x\right)\frac{\partial V_k}{\partial t}=G_k^{\text{axoplasm}}{\cos }^2{\theta }_k\left(z\right)\left[\frac{{\partial }^2{V}_k}{\partial z^2}-{\displaystyle \underset{p=1}{\overset{N}{\sum }}}{W}_{ip}\left(z\right)K_p\left(z\right)\frac{{\partial }^2{V}_p}{\partial z^2}\right]-\cdots$ (15)

$G_{m\mathrm{,}k}{V}_k+I_k^{\text{injected}}+I_k^{\text{fields}}$ (16)

In Equation (15) and (16), $I_k^{\text{fields}}$ is computed keeping in mind that there will be an induced current at the present axon’s node of Ranvier based on the weak impinging field of a distant node and there will be a dependence on the open-close state of the ion channels at the present node. Suppose the state of the node of Ranvier at time t is $l\left(t\right)=\gamma L$ where $0<\gamma <1$, then we assume that, of the available impinging field-induced current, only a fraction ${\gamma }^{\prime }<\gamma$ will be able to enter the node. ${\gamma }^{\prime }$ is strictly less than $\gamma$ since there is a possibility that the currents present at the node oppose the current flow induced by the impinging field. On the other hand there may be a chance that the currents present at the node actually aid the induced current, by being in the same general direction. Thus, as a simple model, if the sign of the impinging current matches that of the nodal current, then ${\gamma }^{\prime }>\gamma$ by a small amount $\epsilon$ and the other way round if the signs are opposite^4,5. Consequently, the current $I_k^{\text{fields}}$ can be written as ${\gamma }^{\prime }\sigma E^{\text{impinging}}$. With these two improvements to the prior formulation, we are in a position to perform computer simulations, which we take up next.

5. Joint Simulation of Tracts and Electric Fields

The two programs developed previously which we utilize for developing the simulations in the present paper are:

1) Ephaptic Program—this program, published in 2019 [20], is able to compute the current-mediated interaction between axons and show the synchronization of action potentials in simple geometries. We refer the reader to that reference.

2) Electric Field Program—this program, developed in 2016-17 [8], is able to compute the electric field generated by an arbitrary number of open ion channels in a node of Ranvier. Using a quasi-static assumption, the generated field was pre-multiplied by a temporal envelope which was trapezoidal in shape, reflecting the profile of current flow through a voltage gated ion channel (the characteristic two picoampere profile [5] ). We refer the reader to [8].

We next summarize how reference [22] was used to merge the two programs described here. Essentially the task is to connect the m, n, $\tilde{h}$, and p variables (also known as the activation and inactivation variables) to the number of ion channels open at a node at any given time. [22] tells us the relation between these variables and the number of open channels with time. To quote:

If, at a given time, the numbers of open $\tilde{h}$, m and n-gates relative to the total numbers of these gates are $\tilde{h}$, m and n, respectively, the expected numbers $N_{Na}$ [... and $N_K$...] of open Na + [.. and K + ...] pores are close to $\tilde{h}{m}^3{N}_{Na}$...

This information is used to pre-multiply the field magnitude generated at a single pore, resulting in the field of the entire node, in its full time-variation, even as the action potential is being computed by the main part of the program. Here, instead of using the quasi-static assumption, the time-variation comes in because $m\left(t\right)$, $n\left(t\right)$, $\tilde{h}\left(t\right)$ and $p\left(t\right)$ are functions of time. Thereby, as the program computes m, n, $\tilde{h}$ and p, it also completes a computation of the sodium field magnitude in a small neighborhood of the node under present consideration. Note that we will actually want the field at various nodal positions of other axons.

5.1. The Program Output

Next, we outline the output of the synthesized program (cf. the programs outlined above). We note that tract-tortuosity, discussed in Section 4, was not considered in this first simulation of the synthesized program; it will be considered in future versions of the program. Instead, we persisted with the formulation discussed in Section 3. Figure 3 summarizes the construction of the program used in this section.

First, we present the modified part of the code, shown in Section 5.1.1. We first discuss the variables used in this algorithm.

1) In this code snippet, bundle1 refers to an entire axon tract, which is commonly also known as an axon bundle.

2) Endlength refers to the length of the axon.

3) Bundle1.N refers to the total number of axons in the bundle.

4) Bundle1.Ranvier refers to the spatial locations of the nodes of Ranvier, within each of the axons in the bundle.

5) Ismember is a method that tests for set inclusion.

6) Activationupdate is a method that simply finds the m, n, h and p values at the next time step, given their present values, the present transmembrane voltage, the temperature-related variable Q, and the time-step size.

Figure 3. Outline of the program used for electric field based interaction simulation.

7) Xlim, plot and drawnow enable the MATLAB functionality of plotting used here.

The first part of the algorithm snippet in the next subsection, that is, the two nested for loops, were a part of the original [8] current-mediated coupling program. These two nested for loops update the m, n, $\tilde{h}$, p variables. As soon as this is done and the loops are exited, we enter into the new part of the program which plots the field at the present time at a pre-selected node of Ranvier on the axon currently under consideration, the j-th axon, via a call to mainporeAug172020.m, a file which deals with the electric field at the ion channels (see the “Nodal Electric Field Computation” block in Figure 3). Thereby, the two programs mentioned at the beginning of this section are interwoven.

Code Output of Overall Program Inclusive of Algorithm Snippet

One of the limitations of this program is that we didn’t use exact interaxonal distances but rather considered the field at a fixed (relatively large) distance from the nodal surface and applied that to the destination node. Considering the exact interaxonal distances would have added computational expense.

• for i = 1:bundle1.endlength;

• for j = 1:bundle1.N;

• if ismember (i, bundle1.Ranvier) % Execute the Ranvier nodal code % Update the concentrations m, n, h, p;

• [axon(j).m(i, t + 1), axon(j).n(i, t + 1), axon(j).h(i, t + 1), axon(j).p(i, t + 1)];

• =activationupdate(axon(j).m(i, t), axon(j).n(i, t), axon(j).h(i, t);

• axon(j).p(i, t), (axon(j).Vcoeff(i, t)), Q, bundle1.delT);

• end % Finish with Ranvier nodes end % Skip myelinated internodes end;

• Figure 1; xlim([1, LTminus1]); plot(1:t, axon(j).m(20, 1:t), “b*”); drawnow;

• Figure 2; plot(1:t, axon(j).n(20, 1:t), “r*”); drawnow;

• Figure 3; plot(1:t, axon(j).h(20, 1:t), “g*”); drawnow;

• mainporeAug172020(axon(j).m(20, t);

• axon(j).n(20, t), axon(j).h(20, t), axon(j).p(20, t), t).

The variable gain pre-multiplies the field strength and is used to control the impact of the field on the other axons. This allows us to model varying conditions such as presence or absence of extracellular organelles, transient changes in local ionic concentrations, and the like. The following Figure 4 shows that when $\kappa =6e-11$ was used for the fields impinging on other nodes, the overall coupling was increased, and led to rapid rise in conduction velocity of the action potentials. The next Figure 5 shows the same conditions as the previous Figure 4, but with a null gain for the field-based coupling. Comparing these two Figure 4 and Figure 5, clearly the network topology, interactions and throughput have changed due to the additional field-based injections. Next, we reduce the inclinations of the axons to 15˚ each, and re-run the high gain setting. We obtain the following Figure 6. Penultimately, we repeat the high gain setting, but this time with ${\theta }_k=5^{\circ }$ for all axons. We obtain the following Figure 7. Finally, we moderately lower $\kappa$ to 5e−11 instead of 6e−11. We perform the simulations for an

Figure 4. Rapid rise in conduction velocity with high $\kappa =6e-11$. N = 3 axons were used, the first one was injected. A moderate current-mediated coupling ratio $r=0.15$ was used. ${\theta }_k={30}^{\circ }$ for all axons k and inter-axonal geometry W as introduced by Chawla, Morgera and Snider was used. At the nodal level, 10 ion channels were used per node, in order to demonstrate the principles on which the program is built, while economizing computational time. The particular sodium channels to be opened at each time step were chosen probabilistically. Only nodes at the same lateral level (facing nodes) interacted across axons via their generated fields.

Figure 5. Null gain for field based interaction. All of the other conditions match the previous Figure 4.

Figure 6. High gain for the field-based coupling and a lower inclination of 15˚. All other conditions remaining the same as in Figure 5. A slight separation between the traces representing the final two nodes is visible on the third axon, as compared with the same axon, same conditions plot of the previous figure. This would indicate a slightly slower conduction velocity on the third axon as compared with the ${\theta }_k={30}^{\circ }$ case.

intermediate ${\theta }_k={10}^{\circ }$ maintaining the same W matrix. We obtain the following Figure 8. We can improve the simulations in future, by allowing every node to influence every node on the remaining axons, with the contribution

Figure 7. High gain setting, but with ${\theta }_k=5^{\circ }$ for the inclinations of the axons. The figure is zoomed in on the transitions from resting to firing. This captures the difference in the three axons’ voltage profiles.

Figure 8. Same as previous figures, but with a moderate ${\theta }_k={10}^{\circ }$ and a lower $\kappa$ of 5e−11. There are clear differences in action potential initiation times at various nodes of the different axons, when compared with Figure 7.

appropriately scaled by the actual distance to the other node. This largely concludes our currents and fields program. In the next section, we will include the effects of the potassium ions and in Section 7 we will look at a new modality, which can perhaps even lead to interaction between tracts themselves, not just within a tract.

6. Including the Potassium Effect

In this section, we show that including the “opposing” effect of potassium is crucial. The effect is opposing because the potassium current is inward whereas the sodium current is outward. Thus, even though the charge on both the ions is of the same sign, the inward flow of the potassium ions is equivalent to a negative outward current and hence the potassium-generated field at the target node will oppose the sodium-generated field. We verify this intuition via simulations. We find that when the potassium “gain” is half of the sodium gain $\kappa$, the output is the same as obtained without considering potassium at all, seen in Section 5. This is shown in Figure 9.

Next, we increase the contribution of potassium current to three-fourths of $\kappa$. The response is shown in Figure 10. As per the figure we see that the response is extended further in time as compared to the 50 percent case. This is significant and shows that it is the potassium current’s field that can lead to a reasonable steady-state wherein the axons are coupled, but yet don’t get overwhelmingly large conduction velocities. This outcome encourages us to further raise the contribution of potassium and re-run the simulation. The result is shown in Figure 11 wherein we have raised the potassium gain to 85 percent of the sodium gain value. Finally, in Figure 12 we raise the gain to 99 percent of the sodium gain and re-run the simulation. We find that, as expected, the point

Figure 9. Three axons are considered, with field-based interaction. The first one is stimulated and the others show a sympathetic response. The potassium electric field’s influence is set at half that of sodium.

Figure 10. Three axons are considered, just as in the previous figure. However, the electric field of potassium has an influence three-fourths that of sodium.

Figure 11. Three field-coupled axons are considered and the first one is stimulated, as in the previous two figures. The potassium gain is set at 85 percent of the sodium gain $\kappa$.

of synchronization⁶ is pushed further towards positive infinity in time. From these studies it is clear that with the potassium current in place we get coupling and, more interestingly, synchronization between the axons, despite the three-dimensional geometry, with the synchronization being more pronounced and early, the weaker the potassium contribution. Since the numbers of

Figure 12. Three axons are considered, with the first one being stimulated, as in the previous figures. Approaching equality, the potassium gain is set at 99 percent that of sodium.

Figure 13. The output of the program for the potassium gain being 95 percent of sodium is shown in this figure. Additional axonal length is appended at the end of the axon. The figure clearly shows a nearly synchronized output after a velocity change. $\kappa$ is 5e−11.

potassium channels at several species’ nodes of Ranvier are smaller than those of sodium channels [9], it is clear that synchronization can take place in axon tracts when they have negligible current-based coupling and full-fledged field-based coupling. With this we conclude our electric field-based coupling studies and enter into a study involving gravitational waves, which are similar to electromagnetic waves in many respects [23].

7. Initial Studies of Gravitational Waves Interacting with the Tract

Here we start with the approach of [24] because it simplifies the interaction problem to a simple metric perturbation. The general problem of derivation of the impact of curvature on the axon tract is left to a future work as it would involve the full machinery of differential geometry [25] on Riemannian manifolds [26]. We encourage our readers to tackle that general derivation. [24] starts with a partial differential equation, the Klein-Fock-Gordon equation, for the purposes of studying Bose-Einstein-Condensates and derives the production of phonons. To study the impact of gravitational waves, the author then states the following invariant spacetime interval:

$\text{d}{s}^2=c^2\text{d}{t}^2-\left(1-h\right)\text{d}{x}^2-\left(1+h\right)\text{d}{y}^2-\text{d}{\tilde{z}}^2$ (17)

which is for a wave with a fixed positive polarization, propagating in the $+\tilde{z}$ direction. For the purposes of our simulations, we will look instead at a wave propagating in the +x direction as that wave is likely to have more interaction with a tract that is oriented in the $+\tilde{z}$ direction. However, when fully tortuous tracts are considered, this choice of direction would become much less important. For such a gravitational wave propagating in the +x direction, the invariant interval becomes, by a cyclic permutation (under the assumption of a spacetime symmetric in the three spatial directions):

$\text{d}{s}^2=c^2\text{d}{t}^2-\text{d}{x}^2-\left(1-h\right)\text{d}{y}^2-\left(1+h\right)\text{d}{\tilde{z}}^2$ (18)

This implies that every partial derivative with respect to $\tilde{z}$ must be replaced with a term that includes a multiplicative prefactor of $\left(1+h\right)$. If the wave is time varying, then h is actually $h\cos \left(2\pi ft\right)$ where f is the frequency of the gravitational wave. Since no other partial derivative actually appears in the GEE equation [20], we limit our analysis to this case. This leads to multiplicative prefactors on terms containing the inclination, within the simulation program.

We can investigate a bit deeper into the background of the gravitational wave setup—beyond [24] —and thus add rigor to the analysis. On a first glance, for linearized gravitational waves (Section 35.5 of [27] ), we need not look at the Einstein field equations (Equation (5.124) in [28] ) or the full equation satisfied by the gravitational wave (Equation (5.153) in [28] ). Following [27], we need only look at the Riemann curvature tensor R which is a quantity like a matrix (Equation (35.10) in [27] ). Based on this tensor, the authors of [27] are able to derive an equation which tells us how the position of one of two particles (particles A and B), placed away from the origin (particle B), will change with time, in terms of the position of the first particle (particle A), placed at the origin. This equation, Equation (35.15) in [27], involves R and the Kronecker delta function.

We may recall that in this analysis the metric g can be expanded into a Minkowski metric term $\eta$ summed with a metric “perturbation” term and higher order terms. The metric perturbation term is related to the gravitational field but not identical with it in general. This is quite convenient for developing an understanding—that these are two separate things. However, in the so-called tranverse-traceless (TT) gauge, which we are using in this analysis, the two turn out to be identical [27].

Coming back to Equation (35.15), if we assume (call this assumption #) that tract dimensions are small such that any distance divided by the speed of light will give a very negligible time duration, then the time-varying strain field, will not depend on location in the tract. Hence we can evaluate Equation (35.15) simply—we will obtain the time varying distance separating the two particles A and B, in terms of the metric perturbation. Taking this as a general principle and applying it to the geometry [20] we have the following. In our case the tract is modeled as a right circular cylinder, increasing in length along the +z axis. Consider an axon in the tract with an inclination $\theta$ to this axis. The coordinate along the axon is s and the projection of s onto the +z axis is z. Thus,

$s\cos \theta =z.$ (19)

If we consider particle A to be at the origin ( $x=y=z=0$ ) and particle B to be along the +z axis, at a distance z from the origin, then we can use the general principle to write down the time variation in z in terms of the metric perturbation. Next, we consider particle B to be at +s along the axon. We would again like to use the general principle to find the time variation in s. If we make further simplifying assumptions:

1) The axon lies in the xz-plane⁷; and

2) The angle $\theta$ changes negligibly with time.

Then we obtain the following description. $s=\left(s_x,s_y,s_z\right)=\left(z\tan \theta ,0,z\right)$, with the consequent length of s being $z\sec \theta$.

Hence, under assumption 2, the time variation in the length of s can be directly obtained from the time variation in z. There will just be a multiplicative pre-factor of $\sec \theta$. But $\sec \theta$ is the inverse of $\cos \theta$ and so under assumptions 1 and 2, the length of s will be in the same relation to z as in Equation (19). And so the time varying length of s will be the time varying length of z, divided by $\cos \theta$ (keeping assumption 2 in mind). The conclusion is that $\text{d}s\left(t\right)=\frac{1}{\cos \theta }\cdot \text{d}z\left(t\right)$. So the same relation holds as before (please see the first relation in Section 2 of [20] ), except that there is time dependence from the metric perturbation.

Based on Equation (35.15), this time dependence is (under assumption #) of the form $z\left(t\right)\sim z\cdot \left(1+h\cos \left(\omega t\right)\right)$. Thus, after a bit more analysis where we study the places where $\cos \theta$ appears in Equation (9) of [20], regrouping and collecting the appearances, we find that wherever in the code we have a $\cos \theta$, it should be multiplied by $\left(1+h\cos \left(\omega t\right)\right)$. With modifications to the program along these lines, we obtain the figures shown in this section. The key assumption is assumption 1—a constraint on the arrangement of axons, which is not too restrictive.

Next we use the thus rederived GEE equation in order to allow x-flowing gravitational waves to have an impact on the tract. In what follows, we present the results of our simulation. Figure 14 shows the state of affairs before any gravitational radiation has arrived. Figure 15 shows what happens when a x-propagating gravitational-wave with a strain value h of 1e−21 is present. Figure 16 shows what is the output if we make the strain level unrealistically high—a value of 0.1—which it might perhaps take near a black hole collision that is emitting the gravitational waves. The conduction velocity changes on axon number 2, but the effect is moderate. From these simulations we can see the effect of piecewise constant strain fields, while gravitational waves actually have time varying strain fields $h\left(t\right)$. Such fields are investigated in Subsection 7.2.

7.1. Discussion on the Time-Frame of the Evolution of Anthropomorphic Brains

We are interested in the evolution [29] of human-like brains [30], traced to its primordial cause at the Big Bang. Thus we take a cosmological perspective [28]. Coming to the Earth stage after the appearance of vertebrates, suppose we could establish that an appearance of a spike in coupled axon number 2 at 2.3 milliseconds was crucial to a particular function of the corresponding nerve tract, such as execution of the “fight-or-flight” response. Then, it is clear that gravitational radiation may play a very important and delicate role in the survival of the vertebrate with a nervous sytem.

Figure 14. The inclination is ${\theta }_k={10}^{\circ }$ for all axons and the same W matrix as in Figure 4. No gravitational radiation is present.

Figure 15. A constant value of h = 1e−21 is used. With such a low strain, typical of gravitational waves observed on the Earth (as indicated by Thorne), there is no change as compared with Figure 14.

Figure 16. Near a source of strong gravitational waves. The time traces are altered. In particular, there are more action potentials visible in the plots, which indicates a rise in conduction velocity. A very large strain h = 0.1 is used.

As suggested by Thorne [31], many of the effects of gravitational waves would be important near the Planck era of the Big Bang and not very significant much after that. Given the above discussed coupling between gravitational waves and action potentials, it may be that consciousness-processing entities such as the brain might not even be possible near the Big Bang (they would be precluded by a host of other reasons too) since the ambient gravitational radiation would scramble neuronal signals to a great extent in that era. Thus the time evolution of the universe from the Big Bang onwards sets a lower limit to the time when conscious entities with stable information-processing nervous systems could have evolved.

It may thus seem that the present paper develops a negative result. However, one of the additional motivations of this paper is to blur the division between what lies inside and outside man. What is inside the brain is deeply tied with what is outside, a realization that would dawn upon a study of this paper, and not be summoned by simply viewing the brain as a black box to which the senses report, with an ego sitting inside which processes that information. Rather the picture painted by this paper is that the ego may be universal in its reach. In a certain sense, Wheeler’s image of the universe as a self-excited “quantum” circuit ( [32] [33] ) is supported by this paper (though we don’t investigate the quantum mechanical aspect) as we propose the “oneness” of the universe, inclusive of the brain of the perceiver, and (hence) “consciousness”.

7.2. Variation of Gravitational Wave Frequency and Its Impact on Relative Time of Initiation of Impulses

In this section, we introduce time variation into the strain and simulate the corresponding modified GEE equation. Figure 17 may be called the geometric frequency response of the tract. If we treat the signal as being encoded in the time

Figure 17. We vary the frequency of a h = (0.1-red, 0.09-blue) strain value gravitational wave and pass it through the tract and record the time gap of sympathetic intiation on axon 3 with respect to axon 2. The presence of the wave scrambles the timing response. Some frequencies are suppressed, and others are magnified. If any important information is time coded, it would be modified. Rate coded information may survive better.

difference of the responses of axons 2 and 3, then we are seeing the response of the system⁸. With certainty, this signal reflects the geometric properties of the tract and encodes its geometric data, such as the axonal inclinations and the W matrix. So it is the geometric signature of the tract. As the figure shows, there exist gravitational wave frequencies for which the signal is magnified and those for which the signal shows abatement.

So in the presence of the metric perturbation of a specific frequency, the tract input (the stimulation on axon 1) leads to a specific time gap between ephaptic coupling on the other two axons, which can be treated as the output. In simpler terms, the sensory or other neural stimualtion to the axon combines with the gravitational wave perturbation signal to generate a specific time lag. Thus local processing as well as gravitational processing are simultaneously encoded in the tract’s response. This may be thought of as cross-modal integration, the two modes being gravitational and electromagnetic, modes of the external world influencing the tract. It is also akin to a transfer function between gravitational radiation and electromagnetism with the tract representing a coupling between these two forces.

In [31], Figure 9.6, the characteristic amplitudes and frequencies of gravitational waves from several postulated periodic sources have been indicated. All the sources studied have very small amplitudes, in the range 10⁻¹⁸ to 10⁻²⁸ and frequencies in the range 10⁻⁶ to 10⁴ Hertz. The frequencies we study are about 10⁻³ Hertz but the amplitudes we need for getting significant effects are very large. It is to be noted however that in works such as [31] or [34], the authors are interested in the minimum amplitude that is detectable at Earth. They are typically not interested in, and do not mention, large amplitudes and their detectability and sources. Thus the absence of mention of amplitudes around 10⁻¹ in the literature, as studied in the present paper, cannot be construed as their physical impossibility of occurrence.

This concludes our study of gravitational radiation and its impact on the tract, and with it, we conclude the physical layer study of the tract. In the next section, we will look at the data or information layer of the tract where it acts as a carrier of information in bits.

8. Channel Capacity of the Tract

The capacity of a neuronal link has been of long-standing interest to neurophysiologists [35]. Ikeda and Manton also have a paper where they look at the capacity of a single neuronal channel [36]. Interesting questions arise when the neuronal link has substructures within it, such as individual axons. Will the inter-axonal connectivities and geometry have an impact? Likely yes, but how strong of an impact on capacity? Within the ambit of our program of further developing the same simulation software as used in the previous sections, we develop a meta-function that is able to numerically find the capacity of an

Figure 18. Variation of tract capacity with axonal inclinations. Three axons are considered with the W matrix as in the previous figures. By changing which ones are stimulated we get 8 possible input conditions. The output is measured at a particular downstream node on each of the three axons. With the inputs and outputs in hand, we can compute the capacity using the Blahut-Arimoto algorithm, by first finding the channel transition matrix via repeated runs of the program, for each of the inclinations studied. Sources of channel randomness are within the tract itself, at the ion channels, as in the previous sections.

ephaptic tract, treating it as a discrete memoryless channel, using the standard Blahut-Arimoto algorithm [37] [38]. We study the variation of the capacity with tract geometry. In particular, as shown in Figure 18, we may vary the angular inclination of the axons, keeping the W matrix constant, and study the resulting tract capacity.

We find an overall increase in tract capacity with the inclinations of the axons. Since the simulations are numerical, this increase is not directly explainable, but may be due to the fact that with increasing inclinations, parts of the axons come much closer to each other and this effect dominates over the other longitudinally opposite parts which tend to move further apart. Perhaps, on coming closer, they are able to allow hopping points for impulses to cross over to their neighbors.

A deeper investigation of this effect would have bearing on our understanding of MS wherein there is a deterioration of information carried by the axon tract with disease progression. Are the inclinations changing in MS? Likely not, but an equivalent transformation must be occurring such as increasing demyelination and the suggestion is that decreasing demyelination is similar in its effect on capacity, to increasing inclinations.

9. Conclusion, Limitations and Future Work

In this work we have viewed the brain as an arena for the interplay of universal forces. In particular, we have looked at axon-axon interaction mediated via electric fields, when the fields influenced the other axons via nodal current injections. We found a dramatic increase in conduction velocity on the coupled axons which was due to the positive feedback in the coupling of the axons. This indicates that electric fields can be a significant mechanism for the change in velocity required for synchronization between axons. With the potassium correction in place, such synchronization was nearly observed.

We also studied what happens when weak and strong gravitational radiation is applied to the axon tract and found effects were present, though small. These latter results simply represent the fact that gravitational interaction is a very weak interaction. Unless the nervous system approaches a singularity, the effects on consciousness via changes in information processing ought to be negligible. Nevertheless, to have a concrete inter-relationship, illumined by simulations, is a step forward in showing the interconnectedness of all systems, particularly, in this case, living systems and the structure of spacetime. Via the gravitational wave, we may even have coupling between different nervous systems. Thus, though the gravitational investigation has yielded largely a negative result, it is nevertheless a pointer to the importance of geometrical studies. Perhaps brain-wide geometrical formulations would yield more interaction with gravitational radiation, even weak waves.

Coming to a discussion of limitations of our work, we may postulate a general stochastic noise related to both the electric field generation mechanism at the ion channels and the gravitational radiation (via stochastic gravity). These together are then biophysical sources of noise in the nervous system—the first one was included here while the second one was not. Another limitation of our work is that none of the simulations presented were carried out for really tortuous tracts. In the presence of tortuosity, even weak gravitational waves may have a significant impact. A further important limitation is that the $m\mathrm{,}n\mathrm{,}\tilde{h}\mathrm{,}p$ variables evolve dynamically with time at a given node of Ranvier. However since in a general relativistic view, spacetime is one unified entity, one ought to include the impact of curvature and gravitational radiation on the evolution of the $m\mathrm{,}n\mathrm{,}\tilde{h}$ and p variables as well. This is left for future investigation. Another study that merits future investigation is to perform the investigations of this paper with axons of varying inclinations rather than identical inclinations.

As a further extension, we can study synchronization between axons in a tract in the presence of gravitational radiation and investigate the entropy amplification between tract input and tract output in comparison to that without gravitational radiation. The difference between the two cases, with and without gravitational radiation, may be denoted as gravitational information, measured in bits. These bits would be part of the information flowing in the nervous system of organisms that have a nervous sytem. They would be the signature of the remotest regions of the cosmos in the information processing of each organism.

In summary, the fields, tortuousity, gravitational waves, and capacity investigations of the present paper have advanced the state of the art in the simulation domain, providing new insights into axon tracts at both the physical layer level as well as the data layer level. These investigations have opened up a number of theoretical questions and inter-relationships among modalities.

Acknowledgements

AC would like to thank Dr. Christopher Passaglia for his review of an earlier draft of this manuscript and Dr. Lav Varshney for suggesting the use of Blahut-Arimoto for channel capacity computation.

NOTES

¹Because the potassium ion is not considered, the results for the impact of the field on conduction velocity (Figures 4-8) show a very sharp rise in the conduction velocity. The additional effect of the potassium field on the propagation is looked at in Section 6.

²In order to avoid repetition, we note here that this new equation is identical with Equations (13) and (14) of the next section, except that we remove the z-dependence in ${\theta }_k\left(z\right)$, $W_{ip}\left(z\right)$ and $K_p\left(z\right)$.

³Please see Figure 1 and associated discussion in [20] for the basic underlying setup which is being assumed in this paper.

⁴In this respect, thinking of the signum function as an activation function, the node itself acts like a mathematical “neuron” (see Figure 7 in Chapter 1.3 of [21] ).

⁵We took the sign variability of ${\gamma }^{\prime }$ into account while running the simulations in Section 5.1.

⁶Here it appears that the synchronization is not sustained and propagation halts after synchronization—but in fact it just reaches the end of the number of available nodes. Please see Figure 13.

⁷The analysis can be repeated similarly for the yz-plane.

⁸A related consideration is: what is the sensitivity of brain information processing to precise timing? Here timing codes would be the more relevant codes (as opposed to rate codes)

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Minsky, M. (1988) Society of Mind. Simon and Schuster, New York.
[2]	Nunez, P.L. and Srinivasan, R. (2006) Electric Fields of the Brain: The Neurophysics of EEG. Oxford University Press, Oxford.
[3]	Landauer, R. (1991) Information Is Physical. Physics Today, 44, 23-29. https://doi.org/10.1063/1.881299
[4]	Chalmers, D.J. (2002) Philosophy of Mind: Classical and Contemporary Readings. Oxford University Press, Oxford.
[5]	Debanne, D., Campanac, E., Bialowas, A., Carlier, E. and Alcaraz, G. (2011) Axon Physiology. Physiological Reviews, 91, 555-602. https://doi.org/10.1152/physrev.00048.2009
[6]	Reutskiy, S., Rossoni, E. and Tirozzi, B. (2003) Conduction in Bundles of Demyelinated Nerve Fibers: Computer Simulation. Biological Cybernetics, 89, 439-448. https://doi.org/10.1007/s00422-003-0430-x
[7]	Pellionisz, A. And Llina, R. (1982) Space-Time Representation in the Brain. The Cerebellum as a Predictive Space-Time Metric Tensor. Neuroscience, 7, 2949-2970. https://doi.org/10.1016/0306-4522(82)90224-X
[8]	Chawla, A. (2017) On Axon-Axon Interaction via Currents and Fields. Ph.D. Thesis, University of South Florida, Tampa.
[9]	Hille, B., et al. (2001) Ion Channels of Excitable Membranes. Vol. 507, Sinauer, Sunderland.
[10]	Selvan, K.T. and Janaswamy, R. (2017) Fraunhofer and Fresnel Distances: Unified Derivation for Aperture Antennas. IEEE Antennas and Propagation Magazine, 59, 12-15. https://doi.org/10.1109/MAP.2017.2706648
[11]	Tønnesen, J., Krishna Inavalli, V.V.G. and Valentin Nägerl, U. (2018) Super-Resolution Imaging of the Extracellular Space in Living Brain Tissue. Cell, 172, 1108-1121. https://doi.org/10.1016/j.cell.2018.02.007
[12]	Horch, K.W. and Dhillon, G.S. (2004) Neuroprosthetics: Theory and Practice. Vol. 2, World Scientific, Hackensack. https://doi.org/10.1142/4987
[13]	Hales, C.G., Grayden, D.B. and Quiney, H. (2011) The Electric Field System of a Macular Ion Channel Plaque. 2011 Annual International Conference of the IEEE Engineering in Medicine and Biology Society, Boston, 30 August-3 September 2011, 294-297. https://doi.org/10.1109/IEMBS.2011.6090077
[14]	Schwinger, J., Deraad Jr., L.L., Milton, K.A., Tsai, W.-Y. and Norton, J. (1998) Classical Electrodynamics. Perseus Books, Reading.
[15]	Plonsey, R. and Heppner, D.B. (1967) Considerations of Quasi-Stationarity in Electrophysiological Systems. The Bulletin of Mathematical Biophysics, 29, 657-664. https://doi.org/10.1007/BF02476917
[16]	Coomar, A., Arntsen, C., Lopata, K.A, Pistinner, S. and Neuhauser, D. (2011) Near-Field: A Finite-Difference Time-Dependent Method for Simulation of Electrodynamics on Small Scales. The Journal of Chemical Physics, 135, Article ID: 084121. https://doi.org/10.1063/1.3626549
[17]	Moritz Schey, H. (2005) Div, Grad, Curl, and All That: An Informal Text on Vector Calculus. WW Norton, New York.
[18]	Polyanin, A.D. and Nazaikinskii, V.E. (2015) Handbook of Linear Partial Differential Equations for Engineers and Scientists. CRC Press, New York. https://doi.org/10.1201/b19056
[19]	Frankenhaeuser, B. and Huxley, A.F. (1964) The Action Potential in the Myelinated Nerve Fibre of Xenopus laevis as Computed on the Basis of Voltage Clamp Data. The Journal of Physiology, 171, 302-315. https://doi.org/10.1113/jphysiol.1964.sp007378
[20]	Chawla, A., Morgera, S.D. and Snider, A.D. (2019) On Axon Interaction and Its Role in Neurological Networks. IEEE/ACM Transactions on Computational Biology and Bioinformatics, 18, 790-796. https://doi.org/10.1109/TCBB.2019.2941689
[21]	Haykin, S. (2010) Neural Networks and Learning Machines, 3/E. Pearson Education India, New Delhi.
[22]	Skaugen, E. and Walløe, L. (1979) Firing Behaviour in a Stochastic Nerve Membrane Model Based upon the Hodgkin-Huxley Equations. Acta Physiologica Scandinavica, 107, 343-363. https://doi.org/10.1111/j.1748-1716.1979.tb06486.x
[23]	Hawking, S. (1966) Properties of Expanding Universes. PhD Thesis, University of Cambridge, Cambridge.
[24]	Schützhold, R. (2018) Interaction of a Bose-Einstein Condensate with a Gravitational Wave. Physical Review D, 98, Article ID: 105019. https://doi.org/10.1103/PhysRevD.98.105019
[25]	Lang, S. (2012) Fundamentals of Differential Geometry. Vol. 191, Springer Science & Business Media, New York.
[26]	Pauli, W. (2013) Theory of Relativity. Courier Corporation, North Chelmsford.
[27]	Thorne, K.S., Misner, C.W. and Archibald Wheeler, J. (2000) Gravitation. W. H. Freeman, New York.
[28]	Weinberg, S. (2008) Cosmology. Oxford University Press, Oxford.
[29]	Darwin, C. (1909) The Origin of Species. PF Collier & Son New York, New York.
[30]	Sarnat, H.G. and Netsky, M.G. (1974) Evolution of the Nervous System. Oxford University Press, Oxford.
[31]	Thorne, K.S., Hawking, S. and Israel, W. (1987) 300 Years of Gravitation. Cambridge University Press, Cambridge.
[32]	Misner, C.W., Thorne, K.S. and Zurek, W.H. (2009) John Wheeler, Relativity, and Quantum Information. Physics Today, 62, 40-46. https://doi.org/10.1063/1.3120895
[33]	Archibald Wheeler, J. (1988) World as System Self-Synthesized by Quantum Networking. In: Agazzi, E., Probability in the Sciences, Springer, Dordrecht, 103-129. https://doi.org/10.1007/978-94-009-3061-2_7
[34]	Nakamura, T., Sasaki, M., Tanaka, T. and Thorne, K.S. (1997) Gravitational Waves from Coalescing Black Hole Macho Binaries. The Astrophysical Journal Letters, 487, Article No. L139. https://doi.org/10.1086/310886
[35]	MacKay, D.M. and McCulloch, W.S. (1952) The Limiting Information Capacity of a Neuronal Link. The Bulletin of Mathematical Biophysics, 14, 127-135. https://doi.org/10.1007/BF02477711
[36]	Ikeda, S. and Manton, J.H. (2009) Capacity of a Single Spiking Neuron Channel. Neural Computation, 21, 1714-1748. https://doi.org/10.1162/neco.2009.05-08-792
[37]	Blahut, R. (1972) Computation of Channel Capacity and Rate-Distortion Functions. IEEE Transactions on Information Theory, 18, 460-473. https://doi.org/10.1109/TIT.1972.1054855
[38]	Joseph, R. (2020) Channel Capacity Using Arimoto-Blahut Algorithm. MATLAB Central File Exchange. https://www.mathworks.com/matlabcentral/fileexchange/32757-channel-capacity-using-arimoto -blahut-algorithm