Introduction

The human brain is a dynamic physical system of unparalleled complexity. While neuroscience has made great strides, many fundamental questions are still unanswered1, including the processes underlying memory formation2, the working principle of anesthesia3, and–most fundamentally–the generation of conscious experience4,5,6. It therefore seems pertinent to explore whether the brain might generate, transmit and store information using other physical modalities than the ones that have been discovered so far.

In the present work we focus on the question whether biophotons could serve as a supplementary information carrier in the brain in addition to the well established electro-chemical signals. Biophotons are the quanta of light spanning the near-UV to near-IR frequency range. They are produced mostly by electronically excited molecular species in a variety of oxidative metabolic processes7,8 in cells. They may play a role in cell to cell communication7,9, and have been observed in many organisms, including humans, and in different parts of the body, including the brain10,11,12,13,14,15. Photons in the brain could serve as ideal candidates for information transfer. They travel tens of millions of times faster than a typical electrical neural signal and are not prone to thermal noise at body temperature owing to their relatively high energies. It is conceivable that evolution might have found a way to utilize these precious high-energy resources for information transfer, even if they were just the by–products of metabolism to begin with. Most of the required molecular machinery seems to exist in living cells such as neurons16. Mitochondrial respiration17,18 or lipid oxidation19 could serve as sources, and centrosomes20 or chromophores in the mitochondria21 could serve as detectors.

However, one crucial element for optical communication is not well established, namely the existence of physical links to connect all of these spatially separated agents in a selective way. The only viable way to achieve targeted optical communication in the dense and (seemingly) disordered brain environment is for the photons to travel in waveguides. Mitochondria and microtubules in neurons have been hypothesized to serve as waveguides22,23,24,25. However, these structures are too small and inhomogeneous to guide light efficiently over significant distances.

Here we propose myelinated axons as potential biophoton waveguides in the brain, and we support this hypothesis with detailed theoretical modeling. These axons are tightly wrapped by a lamellar structure called the myelin sheath, which has a higher refractive index26 than both the inside of the axon and the interstitial fluid outside (see Fig. 1a). This compact sheath could therefore also serve as a waveguide, in addition to increasing the propagation speed of an action potential (via saltatory conduction) based on its insulating property27. There is some indirect experimental evidence for light conduction by axons12,28, including the observation of increased transmission along the axes of the white matter tracts, which consist of myelinated axons29. Myelin is formed in the central nervous system (CNS) by a kind of glia cell called oligodendrocyte. Interestingly, certain glia cells, known as Müller cells, have been shown to guide light in mammalian eyes30,31.

Figure 1
figure 1

3-D schematic representation of a segment of a neuron, and an eigenmode of a cylindrical myelinated axon.

(a) Different parts of a segment of a neuron whose myelinated axon is sliced longitudinally near the end of the segment. The inset depicts the cross section in the transverse plane. Here r and r′ are the inner and outer radii of the myelin sheath, d is the thickness of the myelin sheath, and nmy, nax, and next are the refractive indices of the myelin sheath, the inside of the axon, and the interstitial fluid outside respectively. The compact myelin (shown in red) terminates in the paranodal region near the Node of Ranvier, with each closely apposed layer of myelin ending in a cytoplasm filled loop (shown in light red). (b) Magnitude of the electric field of a cylindrically symmetric eigenmode (λ = 0.612 μm) of a (cylindrical) myelinated axon, with r = 3 μm, and r′ = 5 μm. (c) A vector plot of the electric field showing the azimuthal polarization of the input mode. For clarity in the depiction of the direction of the field at different points, the arrow length is renormalised to the same value everywhere. The adjacent color bar depicts the actual field magnitude. (d,e) Electric field components along the Y (Ey), and Z axes (Ez) respectively.

An interesting feature of photonic communication channels is that they can transmit quantum information as well. The potential role of quantum effects in biological systems is currently being investigated in several areas, including olfaction32,33, avian magnetoreception34,35, and photosynthesis36,37. There is also growing speculation about the role of fundamental quantum features such as superposition and entanglement in certain higher level brain functions15,24,38,39,40. Of particular relevance is the “binding problem” of consciousness, which questions how a single integrated experience arises from the activities of individual molecules in billions of neurons. The answer to this question might be provided by quantum entanglement41, where the whole is more than the sum of its parts in a well-defined physical and mathematical sense.

The main challenge in envisioning a “quantum brain” is environmental decoherence, which destroys quantum effects very rapidly at room temperature for most physical degrees of freedom42. However, nuclear spins can have coherence times of tens of milliseconds in the brain43,44, and much longer times are imaginable40. Long-lived nuclear spin entanglement has also been demonstrated in other condensed-matter systems at room temperature45. A recent proposal on “quantum cognition”40 is based on nuclear spins, but relies on the physical transport of molecules to carry quantum information, which is very slow. In contrast, photons are well suited for transmitting quantum information over long distances, which is why currently envisioned man-made quantum networks rely on optical communication channels (typically optical fibers) between spins46,47.

Efficient light guidance therefore seems necessary for both classical and quantum optical networks in the brain. Is this possible in myelinated axons with all their “imperfections” from a waveguide perspective? In an attempt to answer this question, we have developed a detailed theoretical model of light guidance in axons. We show in the next section that the answer seems to be in the affirmative.

Results

Introduction to our approach

We use Lumerical’s software packages FDTD (Finite Difference Time Domain) Solutions and MODE Solutions for numerically solving the three dimensional electromagnetic field equations in various scenarios to elucidate the waveguiding characteristics of myelinated axons. For the majority of our simulations, we take the refractive indices of the compact myelin sheath, the axon and the fluid outside to be 1.44, 1.38 and 1.34 respectively (see Fig. 1a), consistent with their observed values26,48,49. This index contrast allows guided modes of light inside the myelin sheath. Although there are many scatterers both inside the axon (cell organelles, e.g. mitochondria, and endoplasmic reticulum, lipid vesicles, and filamentous structures, e.g. microtubules, and neurofilaments) and outside (e.g. microglia, and astrocytes), modes confined primarily in the myelin sheath will effectively not see these structures. The modes are primarily confined in a waveguide if its dimension is close to or larger than the wavelength of the light. Myelinated axons in the brain greatly differ in length and caliber. The short axons of the interneurons are only ~1 mm long, while the longest axons can run through almost the whole length of the brain with numerous branches. Their diameters range from 0.2 microns to close to 10 microns50. We shall assume the g-ratio (the ratio of the radius of the axon, r and the outer radius of the myelin sheath, r′) equal to 0.6 for the majority of our simulations, close to the experimental average51.

Biophotons have been observed in the wavelength range 0.2 μm–1.3 μm. Most proteins (including the proteins in the myelin sheath) strongly absorb at wavelengths close to 0.3 μm. To avoid absorption, we shall consider the transmission of light with wavelengths above 0.4 μm. For different axon calibers, we send in light at different wavelengths, ranging from 0.4 μm to the thickness of the myelin sheath (denoted by d), or 1.3 μm (the upper bound of the observed biophoton wavelength), whichever is smaller. This ensures good confinement in the myelin sheath to limit interactions with the inhomogeneous medium inside and outside the axon (see Supplementary Information). We call this upper wavelength bound the longest permissible wavelength (λmax). The shortest permissible wavelength (λmin) for all simulations is 0.4 μm. In addition to λmax, and λmin, we choose an intermediate wavelength corresponding to the central permissible frequency (mid-frequency of the permissible frequency range), denoted by λint. In a single simulation, FDTD calculates the input mode at λint, and sends light at different wavelengths with the same spatial mode profile (see Methods). Note that for the thinnest axons considered, λmax = λint = λmin = 0.4 μm (d = 0.4 μm, too, for good confinement).

Guided modes in myelinated axons

Let us pick the thickest axon in our analysis (r = 3 μm, and r′ = 5 μm) to explain a few relevant details associated with the guided modes. Figure 1b shows the magnitude of the electric field of a cylindrically symmetric eigenmode of this structure for the wavelength 0.612 μm, obtained using the Finite Difference Eigenmode (FDE) solver in Lumerical’s FDTD Solutions. The electric field is azimuthally polarized (see Fig. 1c–e), similar to the TE01 mode of a conventional fiber52, where the refractive index of the core is higher than that of the cladding. Azimuthal polarization would prevent modal dispersion in the birefringent myelin sheath, whose optic axes point in the radial direction53. There are hundreds of other guided modes for this thickness of myelin sheath. Photons generated by a realistic source in the axons could couple to all these modes, with various coupling coefficients. However, for the sake of simplicity (and lack of knowledge of realistic photon emission characteristics by potential sources), we pick a single mode and study its transmission.

Next, we discuss transmission in the presence of several optical imperfections (e.g. discontinuities, bends and varying cross-sections). The theory of various imperfections in optical fibers is developed with long distance communication in mind, which requires very small imperfections, and focuses on the conventional fiber geometry. Since the myelin sheath based waveguide does not pertain to either of these conditions, we shall mostly deal with explicit examples. We simulate short axonal segments as the computational resource requirements for FDTD are very high, and extrapolate the results for the full length of an axon.

Nodal and paranodal regions

The myelin sheath is interrupted at almost regular intervals by the ‘Nodes of Ranvier’, leaving the axon bare for approximately 1 μm54. The lamellae, whose fusion and apposition leads to the formation of the compact myelin, terminate near the nodes in the paranodal region such that each lamella ends in a loop filled with dense cytoplasm (see Fig. 1a). Many of these cytoplasmic loops are attached to the axonal membrane. For a thin myelin sheath, the paradonal region is almost ordered, with the innermost lamella terminating first, and so on, but for thicker sheaths, many cytoplasmic loops do not reach the axonal surface, but terminate on other loops. The length of the paranodal region should, however, depend roughly on the the number of myelin lamellae, and hence on the thickness of the myelin sheath. We shall call the ratio lparanode/d, the p-ratio, where lparanode is the length of the paranode in one internode (the axon between two consecutive nodes) and d is the thickness of the myelin sheath as defined earlier; p-ratios around 5 are realistic55.

Figure 2a shows our model for the region of the axon close to the node (see Methods section for detailed discussion on the model), and Fig. 2b depicts the magnitude of the electric field in the longitudinal direction (along the length of the axon), as a cylindrically symmetric input mode crosses this region. We call this EFPL (Electric Field Profile in the Longitudinal direction). We vary the length of the paranodal region for an axon with r′ = 5 μm in Fig. 2c, and observe the corresponding modal transmission (power transmission in the guided modes) up to a wavelength away from myelin sheath boundaries (see Methods and Supplementary Information) for 3 different wavelengths, 0.40 μm, 0.61 μm, and 1.30 μm. We interpret the results in terms of beam divergence and scattering, which are the two primary sources of loss here. Shorter wavelengths diverge less, but scatter more. For a short paranodal region (p-ratio = 2.5), shorter wavelengths have a higher transmission, but for longer paranodal regions, longer wavelengths fare better because scattering becomes a more powerful agent of loss than divergence as the length is increased. In general, the transmission drops for all the wavelengths as the p-ratio increases, although the trend is less clear for the longest permissible wavelength.

Figure 2
figure 2

Nodal and paranodal regions.

(a) Longitudinal cross-section of our 3-D model of the nodal and the paranodal regions. For this example, the radius of the axon including the myelin sheath, r′ = 5 μm, g-ratio = 0.6, and the length of the paranode, lparanode = 5 μm. (b) Magnitude of the electric field in the longitudinal direction (EFPL) as a cylindrically symmetric input mode with wavelength 0.612 μm crosses the region. (c) Transmission for an axon with r′ = 5 μm, as a function of the p-ratio, defined as lparanode/d. (df) Transmission as a function of the axon caliber for different wavelengths and different paranodal lengths. The number of cytoplasmic loops, and the microtubules in the paranodal region are kept in proportion to the thickness of the myelin sheath, and the volume of the paranodal region respectively (see Methods).

In Fig. 2d–f, we simulate the nodal region for 5 different axon calibers, several different wavelengths, and different pratios. In general, the greater the mode volume, the less is the divergence. So, a mode with a larger mode volume (corresponds to thicker myelin sheath) should diverge less for the same paranodal length. Here, however, we are dealing with ratios (lparanode/d), rather than absolute values of the lengths, making intuition slightly difficult. Still, in Fig. 2d, we see that the most loosely confined modes (λmax) crudely follow this intuition, and transmission increases for thicker axons. For a fixed axon caliber, the transmission does not depend on the paranodal length in a well-defined way. One possible explanation for this feature is the unconventional nature of the waveguide itself. The long wavelengths mostly suffer loss because of divergence. However, in these waveguides, not all the light that diverges is lost. There is a possibility of a fraction of the light diverging into the axon to come back into the myelin sheath at the end of the paranodal region. This sometimes even increases the transmission in the myelin sheath for longer paranodal regions.

In Fig. 2e,f, for p-ratio = 2.5, the trend follows the intuition based on divergence. Increase in myelin thickness leads to better confinement, and less divergence. However, for larger pratios, the trend almost reverses, and thicker axons perform worse than thin ones. In these cases, scattering becomes more relevant than divergence, and longer paranodal regions lead to greater scattering.

Even with such a sudden discontinuity in the sheath, we find that transmission can still be fairly high. To summarize, if the p-ratio is small (~2.5), well confined modes (shorter wavelengths) yield higher transmission, whereas loosely confined modes fare better for larger p-ratios. Thicker axons are usually better than the thinner ones for smaller pratios (~2.5) at all wavelengths. However, for shorter wavelengths and larger pratios (~5 or greater), thinner axons have higher transmission. We verified that the transmission after multiple paranodal regions can be approximately predicted by exponentiating the transmission through one (see Supplementary Information).

Bends

Optical power from the eigenmodes of a straight waveguide leaks out on encountering bends. Bends of constant curvature have eigenmodes which can propagate with minimal loss, but axons typically change their curvature along their length. These bent–modes (eigenmodes for circular bends) are more lossy than the straight–modes (eigenmodes for straight structures) for changing curvature. Therefore, an appropriate way to quantify the bend losses for an arbitrary axon path will be to incident the straight–mode in a waveguide with continuously varying curvature, and observe the transmission in the myelin sheath at the other end. We choose a sinusoidal waveguide since it has alternating regions of positive and negative curvatures, and can thus serve as a prototype for any arbitrary contour. Figure 3a is an example for an axon with radius 0.6 μm, and Fig. 3b shows the EFPL as a straight–mode passes through.

Figure 3
figure 3

Bends.

(a) The geometry of a sinusoidally bent waveguide. For this example, r′ = 1, and the amplitude (A) and wavelength of the cosine function (l) are 5 μm and 100 μm repectively. (b) EFPL as the input mode with wavelength 0.4 μm crosses the region. (c) Transmission as a function of the change in curvature, Δκ, for different wavelengths in an axon with r′ = 5 μm (Δκ is varied by varying A). (df) Transmission as a function of the axon caliber for different wavelengths and different Δκ.

Bending losses for conventional S bend waveguides (half a cosine function) depend most strongly on the change of curvature56. We therefore plot total power transmission (calculated by integrating the real part of the Poynting vector of the output light directly across the required area, and dividing it by the source power) up to a wavelength away from the myelin sheath boundaries (see Methods and Supplementary Information) as a function of the change of curvature, Δκ = 4Ak2 (k is the wavenumber of the sinusoidal function) for 3 different wavelengths in Fig. 3c (r′ = 5 μm). A shorter wavelength is better confined and therefore yields higher transmission. Fig. 3d–f compare transmission for axons of different calibers. Note that we calculate Δκ of the curve passing through the central axis of the axon. But the inner part of a bent axon has a larger curvature than the outer part at each point. Such a difference becomes particularly important for thicker axons, since they see a larger effective change of curvature than thinner axons, and suffer more loss for the same Δκ. For Δκ ~0.024 μm−1, almost all the permissible wavelengths are guided with negligible loss for all axon calibers discussed. We assume that in a typical axon, regions of large curvature do not exist for considerable length (which seems justified57) and Δκ is a good parameter to quantify the bend inhomogeneity. Some of the axons (r′ ~1 μm) of length ~1 mm in the images of ref. 57 appear relatively straight with Δκ < 0.05 μm−1, which would yield greater than 90% transmission.

Varying cross-sectional area

The thickness of the myelin sheath is not uniform all along the length of the axon. We vary d randomly according to a normal distribution, assuming the correlation length in the roughness of the myelin sheath boundary to be 5 μm–10 μm (see Methods). The mean of the distribution is in close agreement with that predicted by the g-ratio, and the standard deviation (s.d.) of the distribution is varied. Figure 4a shows the longitudinal cross-section for one such simulation (r = 1.8 μm, length of the axonal segment is 50 μm, and the s.d. is 30% of the average thickness of the myelin sheath). Figure 4b shows the EFPL for input light with λ = 0.48 μm. In Fig. 4c, we observe that, in general, a more random distribution of the radius suffers a greater loss (for all wavelengths), and shorter wavelengths transmit slightly better. Figure 4d–f compare the total power transmission (up to a wavelength from the myelin sheath boundaries) in axons of different calibers. Thinner axons can tolerate greater percentage-inhomogeneity, suggesting a closer dependence on the absolute value of the inhomogeneity. All the axons have close to unit efficiency in transmission for less than 10% variation in radius. Extrapolation for transmission in a longer segment of the axon is straightforward. One can exponentiate the transmission fraction by the number of 50 μm segments in the axon. Longer correlation lengths would yield better transmission for the same s.d, while significantly shorter correlation lengths would strongly scatter the mode. Some of the axonal segments (length ~5 μm) of thin axons (r ~ 1 μm) are within this inhomogeneity, as seen in the images of ref. 58. We did not find suitable images of thicker myelinated axons, and longer segments from which a more realistic estimate of this particular inhomogeneity could be extracted.

Figure 4
figure 4

Varying cross-sectional area.

(a) The geometry of a myelinated axon where the cross sectional area of the myelin sheath varies smoothly along the longitudinal direction. For this example, the mean radius of the axon with the myelin sheath is 3 μm and the standard deviation (s.d.) of the variation of the myelin sheath’s radius is 0.36 μm. (b) EFPL as the input mode with wavelength 0.48 μm crosses the region. (c) Transmission as a function of the s.d. of the variation in the myelin sheath’s radius for different wavelengths (the mean radius of the axon with the myelin sheath is 5 μm). (df) Transmission as a function of the mean radius of the axon including the myelin sheath for different wavelengths and different s.d. of the variation of myelin sheath’s radius.

Non-circular cross section

Axons can have quite arbitrary transverse cross-sectional shapes, and the ensheathing myelin partly imbibes that shape58. We give an example in Fig. 5a and the corresponding EFPL when an eigenmode for a circular cross-section (λ = 0.612 μm, r = 3 μm, and r′ = 5 μm) is incident on it in Fig. 5b. In this example, the points along the cross-sectional circumference of the axon are generated randomly according to a normal distribution with a mean value 3 μm and a standard deviation 0.4 μm (13.33% of the axon radius). The myelin sheath is an approximate parallel curve drawn at a perpendicular distance of 2 μm (so that the average g-ratio = 0.6) surrounding the axon.

Figure 5
figure 5

Non-circular cross section of the axon and myelin sheath.

(a) An example of the cross-section of a myelinated axon. The mean distance of the points along the circumference of the axon from its center is 3 μm and the s.d. is 0.4 μm. The outer boundary of the myelin sheath is a parallel curve drawn at an approximate separation of 2 μm from the axonal boundary. (b) EFPL as a cylindrically symmetric eigenmode for a circular cross-section (λ = 0.612 μm, r = 3 μm, and r′ = 5 μm) passes a straight waveguide with this non-circular cross-section. (c) Transmission as a function of the s.d. of the distance between the points on the circumference of the axon and a circle of radius 3 μm for different wavelengths.

Figure 5c shows the total power transmission (within a wavelength of the myelin sheath boundaries) in the myelin sheath for different shapes and different wavelengths in a 100 μm long structure. As expected, transmission drops for all wavelengths as the cross-section becomes more random. Images in ref. 58 show many axons with less than 10% inhomogeneity in the cross-sectional shape. If the axon and myelin sheath do not change the cross-sectional shape substantially along their length, there will be almost no more loss, as coupling loss is the primary source of loss here (rather than the propagation loss). However, if the cross-section changes significantly, there will be propagation loss as well. We model this by twisting an axon, such that it starts out with an elliptical cross-section, interchanges the minor and major axes midway of the simulated segment (25 μm), and reverts to its original shape at the end of the segment (see Supplementary Figure 2).

Other imperfections

In addition to the sources of loss discussed above, there can be several other imperfections, the most significant of which is the cross-talk between axons. Light in a myelinated axon would not leak out significantly, even if placed in direct contact with cells of lower refractive indices than the myelin sheath. However, if two or more myelinated axons are placed very close to each other (side by side), then light leaks out from one to the other (see Supplementary Information). Cross-talk can be interpreted both as a loss and a coupling mechanism between axons in a nerve fiber. In general, the axons should be a wavelength apart to prevent cross-talk, which seems realistic from some of the images in ref. 58.

The other imperfections that we considered do not affect the transmission significantly. The refractive indices of the axon, the myelin sheath, and the outside medium, were taken to be constants for the simulations so far. Next, we varied the refractive indices of the axon, and the myelin sheath, both transversely and longitudinally (with a correlation length for the random longitudinal variation of the refractive index ~5 μm–10 μm), keeping the mean the same as the one used so far, and a s.d. of 0.02 (typical variation as expected from refs 26 and 48) for a few of the simulations. We observed no significant changes in the transmission (typically less than 1%). Moreover, there can be astrocytes and other glia cells in the nodal region close to the axon. As light crosses this region from one internode to the other, it will pass through these cells. We modeled them as spheres with radii varying from 0.1 μm to 0.3 μm, and refractive index 1.4, filling up one third of the volume of the nodal region outside the axon (expected from the images in ref. 58). The transmission through a node of Ranvier increased slightly (~2%) for the thinnest axons, while it stayed almost unchanged for the thickest ones. It is worth noting that in some regions of the internode (apart from the paranodes), the myelin sheath might be unusually inhomogeneous (for instance, by decompaction of the lamellae) leading to additional scattering losses.

Absorption

In biological tissues, and more so in the brain, scattering of light, rather than absorption, is the main source of attenuation of optical signals59. To our knowledge, the absorption coefficient of the myelin sheath has not been measured experimentally. We can only infer it indirectly with limited accuracy. The average absorption coefficient in the white matter decreases almost monotonically from ~0.3 mm−1 to ~0.07 mm−1 for wavelengths 0.4 μm to 1.1 μm60. But myelin can not be responsible for the majority of the absorption since grey matter (almost devoid of myelin) has comparable absorption coefficients60. It is likely that light sensitive structures (e.g. chromophores in the mitochondria) are the main contributors to the absorption. Another way to infer myelin’s absorption coefficient is to look at the absorption of its constituents, i.e. lipids, proteins and water. Mammalian fat shows an absorption coefficient less than 0.01 mm−1 for the biophotonic wavelength range61. Water has similar absorption coefficients. Most proteins have a strong resonance peak close to 0.28 μm with almost negligible absorption above 0.34 μm, and the proteins in the myelin (e.g. myelin protolipid protein, and myelin basic protein) behave similarly62. Thus, absorption in myelin for the biophotonic wavelengths seems negligible (over a length scale of ~1 cm), based on the data of its constituents. Only a direct measurement could tell us more.

Attainable transmission

We discuss a few examples to estimate the attainable overall transmission. The internodal length is typically equal to 100–150 times the axonal diameter51,63. For an axon with r = 3 μm (r′ = 5 μm), internodal length = 1 mm, wavelength of input light = 1.3 μm, s.d. for varying area = 2.5%, Δκ = 0.039 μm−1, s.d. for non-circularity in cross-section shape = 13.33%, separation from the nearby axons = 1 μm, and p-ratio = 7.5, the transmission after 1 cm would be ~31% (see the Methods for the procedure). However, if the wavelength of input light = 0.61 μm, p-ratio = 2.5, and all the other parameters are kept the same, the transmission could be ~82%. A thinner axon with r = 1.8 μm (r′ = 3 μm), internodal length = 500 μm, wavelength of light = 1.2 μm, p-ratio = 7.5, s.d. for varying area = 20%, separation from other axons = 1.2 μm, and Δκ = 0.039 μm−1 would yield ~3% transmission after 1 cm. However, there are neurons in the brain whose axons are ~1 mm long64 (e.g. the local interneurons). If we take a 2 mm long axon, then the transmission for the 3 examples discussed above would be ~78%, ~96%, and ~46% respectively. The predominant loss for these examples is in the paranodal regions. Sources and receivers would need to be located close to the ends of the myelinated sections of the axon to reduce coupling losses. Let us note that photons could travel either way (from the axon terminal up to the axon hillock or the other way round) in an axon.

Attainable communication rates

One potential challenge for the use of biophotons for inter-neuron communication is the fact that biophoton emission rates per neuron seem to be quite low. In ref. 12, the authors count the number of biophotons emitted per minute by a slice of mouse brain, using a photodetector placed at a certain distance away, after exciting the neurons with glutamate, the most common excitatory neurotransmitter. Substituting the relevant experimental parameters, the estimated rate of biophoton emission is about 1 photon per neuron per minute, which is 1–2 orders of magnitude lower than the average rate of electrochemical spikes65. However, this estimate has significant uncertainty. On the one hand, the brain slice is strongly stimulated by glutamate, so the estimate might be high. On the other hand, only the scattered photons are counted. If there are photonic waveguides in the brain, most of the photons propagating in these waveguides would likely be absorbed in the brain itself rather than being scattered outside, so the estimate could be much too low. It should also be noted that the emission rates could be very different depending on the specific neuron or neuron type. Taking the above estimate at face value for the sake of the argument, such low photon rates could still be relevant. Given that there are about 1011 neurons in a human brain, there would still be over a billion photon emission events per second. This could be sufficient to transmit a large number of bits, or to distribute a large amount of quantum entanglement. In this context, it is worth keeping in mind that psychophysical experiments suggest that the bandwidth of conscious experience is less than 100 bits per second66,67. From a quantum perspective, it is known that the behavior of even moderate numbers of qubits (e.g. one hundred) is impossible to simulate efficiently with classical computers68.

Proposals to test the hypothesis

There is some indirect evidence of light guidance in axons12,28,29. A relatively simple way to look for evidence of guidance in-vitro might be to illuminate a brain slice (preferably with many suitable myelinated axons) from one end, and observe for bright spots corresponding to the open ends of suitable myelinated axons at the other end. Optimal illumination might be required to overcome the strong scattering and absorbance by the surrounding medium. For a more conclusive test, one could isolate a neuron with the necessary thickness of the myelin sheath, and small enough inhomogeneity, suspend it in a suitable solution to keep the cell alive for some time, and try to couple one of the guided modes into the axon, similar in spirit to the procedure adopted for verifying light guidance by Müller cells30. We could couple the light in close to the axon terminals, as real sources are suspected to be present there12. To inject a guided mode in the myelin sheath and verify its guidance, one might need to decapitate the axon near the terminal and hillock regions, couple the mode directly in the myelin sheath, and observe the intensity (and if possible, the modal structure) of light emanating from the other end quickly, since the cellular properties start to change soon after death. Evanescent coupling and readout of light is another option.

For an in-vivo test of light guidance, one might first try to prove the presence of photons in the myelin sheath. One could inject a light sensitive chemical (e.g. AgNO3) either in the cytoplasmic loops in the paranodal region directly or in the oligodendrocytes, which would then circulate the chemical in the cytoplasmic loops, and possibly some to the myelin too. Light will activate the oxidation of Ag+ to Ag, which should be visible as dark insoluble granules. This is similar in spirit to the development of photographic films, and the in-situ biophoton autography (IBA) technique28.

Another interesting type of in-vivo tests would involve the artificial introduction of sources and detectors into living neurons. Fluorescent molecules or nano-particles could serve as sources, and also as detectors, if their fluorescence can be triggered by the absorption of photons from the molecule or nano-particle that serves as the source69. An alternative possibility for the insertion of detectors may be provided by optogenetics70, where specific kinds of neurons are genetically modified to produce proteins which can function as light sensitive ion-channels (e.g. channel rhodopsin). If we could embed these proteins specifically in the axonal membrane near a terminal end of the myelin sheath, or in the membranes of the cytoplasmic loops in the paranodal region at an end, and detect photons produced by an artificial source at the other end of the axon, we could verify the light guidance hypothesis. It is interesting to note that there is an increase in oligodendrogenesis and myelin sheath thickness near these genetically modified neurons when stimulated by light71. Do the axons adapt themselves for better light guidance too (in addition to electrical guidance) by adding sufficient layers of myelin?

The final type of test would involve identifying naturally occurring sources and detectors in real neurons, and showing that photons are guided from the sources to the detectors. To our knowledge, photon emission has not yet been studied at the level of individual neurons. Photon measurements have been done macroscopically, counting only the scattered photons10,11,12,13 (neglecting those which are guided or absorbed). It would be important to precisely pinpoint the sources of these photons and to characterize their wavelength and emission rates. This may be possible by enhancing the emission rates through nanoantennas72. It would also be very interesting to study the photon detection capabilities of potential natural detectors, such as centrosomes20 and chromophores in mitochondria21, ideally at the single-photon level. There may be other potential detectors that are yet to be discovered, e.g. light-sensitive proteins similar to channel rhodopsin used in optogenetics70.

Discussion

We have shown that light conduction in a myelinated axon is possible even with realistic imperfections. We have proposed experiments to verify the key aspects of our hypothesis. We now briefly mention several related fundamental questions.

If photons are to serve as quantum communication links between nuclear spins, one also needs to explain how the photons and spins would interface with each other. Researchers in spin chemistry73 have discovered various ways in which electron and even nuclear spins can influence chemical reactions, which can also involve photons. A well-known biological example is provided by cryptochrome proteins, which can be activated by light to produce a pair of radicals with correlated electronic spins, which are suspected to be involved in bird magnetoreception (the ability to perceive magnetic fields)34. Recent theoretical work suggests that interactions between electron and nuclear spins in cryptochromes are important for explaining the precision of the magnetoreception35. Cryptochromes are found in the eyes of mammals too (including humans), and they are also magnetosensitive at the molecular level74. Similar proteins, if present in the inner brain regions, might act as an interface between biophotons and nuclear spins.

In order to connect individual quantum communication links to form a larger quantum network (allowing for the creation of entanglement between many distant spins), the nuclear spins interfacing with different axons would have to interact coherently, which might require close contact. The existence of synaptic junctions between individual axons is particularly interesting in this context.

Concerning the potential relevance of (classical or quantum) optical communication between neurons for consciousness and the binding problem, an interesting anatomical question would be whether brain regions that have been implicated in consciousness75, such as the claustrum76,77, the thalamus, hypothalamus and amygdala78, or a recently identified “hot zone” in the posterior cerebral cortex75 have myelinated axons with sufficient diameter to allow light guidance.

From a medical perspective, an active role of the myelin sheath as an optical waveguide, in addition to the conventional role as an insulating layer, might also enable us to understand the causes of the diseases associated with it (e.g. multiple sclerosis79) better, and help conceive and design subsequent treatments.

Given the potential advantages of optical communication, one might wonder why evolution did not shift exclusively to this modality. Did evolution reach a local optimum of some kind?

If optical communication along myelinated axons is indeed a reality, this would reveal a whole new aspect of the brain, with potential impacts on many fundamental questions in neuroscience.

Methods

Software packages

We use Lumerical’s FDTD Solutions, and Lumerical’s MODE Solutions for all our simulations. Both these software packages use the Finite Difference Eigenmode (FDE) solver to generate the propagation modes for different waveguide geometries. FDE solves Maxwell’s equations for the eigenmodes on a cross-sectional mesh using the finite difference algorithm80. Finite Difference Time Domain (FDTD) method solves Maxwell’s equations in time-domain on a discrete spatio-temporal grid formed by Yee cells81. Since FDTD is a time domain technique, it can cover a wide-frequency range in a single simulation. We use this feature to study the dependence of light guidance on the wavelength of the input light. But one has to be careful in the interpretation of the results. The two main areas of concern are the meshing accuracy of the simulation, and the change of the beam profile with wavelength. Lumerical’s meshing algorithm refines the mesh for smaller wavelengths while leaves it coarse for the larger ones. We manually increase the mesh accuracy for all our large wavelength simulations to keep the number of Yee nodes almost constant for different wavelengths. We have some tiny structures in our simulations (e.g. the microtubules), which need to be included in the mesh. We ensure that they are included by increasing the mesh accuracy to a setting such that the results converge for finer mesh. The variation of the beam profile with wavelength requires careful analysis too. When we select a wide wavelength range, e.g. 400 nm to 1300 nm (equivalently 750 THz to 231 THz), and calculate the eigenmodes of the structure, FDTD calculates the eigenmode at the wavelength corresponding to the central frequency, which is 612 nm (equivalently 490 THz) for this example. It injects light at different wavelengths but with the same spatial field profile. However, the mode-profiles for different wavelengths can differ substantially. Different kinds of waveguide imperfections need different analysis methods to account for this error, and we shall address this point individually for each one. We also ran multiple simulations (narrow wavelength sources at different wavelengths), where we send in the exact eigenmodes, and ensure that the results converge with that for a single simulation and a wide wavelength range.

Paranodal region

The paranodal region is modelled by carving part of a paraboloid out of the cylinder comprising the axon and the myelin sheath, closely imitating their real geometry55. This part of the paraboloid is generated by revolving a segment of a parabola about the axis parallel to the length of the axon. This segment starts at the end of the paranodal region away from the node and terminates at the node. The general equation of the segment is , where r, d, lparanode, x, and y are the radius of the axon, the thickness of the myelin sheath, the paranodal length, the coordinate along the axis of the axon, and the coordinate perpendicular to the axis respectively. The paranodal region is divided into many cytoplasmic loops, modeled by the compartments between concentric rings of increasing radius as one approaches the node. This is in accordance with the fact that the lamellae close to the axonal membrane terminates first and the most distant lamellae terminates last. The thickness of a ring is 10 nm, which is the typical thickness of the cell membrane. The number of these loops equals the number of the lamella in the compact myelin (average separation between adjacent lamellae is 20 nm82). The microtubules in the cytoplasmic loops are generated randomly according to a uniform distribution, and placed transverse to the axon axis. They are concentric cylinders with inner and outer diameter equal to 13 nm and 6 nm respectively. The number of the microtubules is proportional to the volume of the paranodal region. The volume fraction of the microtubules (with respect to the paranodal region) is kept at 1.33% which is a typical value of the volume fraction inside the axon83. The refractive indices of the cytoplasmic loops, the cell membrane, and the microtubules are taken to be 1.38, 1.50, and 1.50 respectively, close to their expected values84,85.

We obtain the electric field profile after a paranodal region, and expand it in an eigenbasis of the guided modes. Each time, we truncate this basis manually (for different axon calibers, and different wavelengths) in Lumerical’s MODE Solutions, neglecting higher order modes (with effective refractive indices, neff < 1.34), almost all of which are lossy. Thus, our basis comprises of guided modes with neff between 1.44 and 1.34. The input mode is antisymmetric about two orthogonal axes in the cross-sectional plane (see Fig. 1d,e). We label these axes Y and Z, with the origin at the center of the axon. Then Ey is antisymmetric about the Y axis, and Ez is antisymmetric about the Z axis, where Ey, and Ez are the real parts of the Y and Z components of the electric field E respectively. Since the input mode is antisymmetric, and the structure is cylindrically symmetric, the guided modes will all be antisymmetric about the central axes. We therefore work in the subspace of antisymmetric modes, and expand the electric field profile in the basis of this subspace. Our waveguide permits a few guided modes primarily confined inside the axon. But the medium inside the axon (and outside the myelin) is expected to be scattering (see Supplementary Information), and we neglect the fraction of power residing in regions beyond a wavelength of the myelin sheath boundaries. We just integrate the real part of the Poynting vector (with the electromagnetic fields corresponding to the guided portion of the output light) across the myelin sheath up to a wavelength from its boundaries (see Supplementary Information for the mathematical expressions). This is an approximate way of expanding the electric field profile in a basis of modes confined strongly in the myelin sheath. We choose to include the power within a wavelength of the myelin sheath boundaries to account for the evanescent fields and a few guided modes which are still very close to the myelin sheath, and are not strongly affected by the inhomogeneities inside and outside.

To account for the change in mode profile with wavelength, we expand the input mode (calculated at the central permissible frequency) in the basis of the guided modes at a particular wavelength. For shorter wavelengths, almost the entire power resides in a superposition of the guided modes (typically greater than 99.5%), but for longer wavelengths, the percentage of power in the guided modes can be significantly lower (e.g. for d = 2 μm, the expansion of an input mode in a basis of the modes at 1.3 μm yield 97.11% coupling). So, we divide the output power (after the paranodal region) for longer wavelengths (obtained by integrating the real part of the Poynting vector with the electromagnetic fields corresponding to the guided portion of the output light across the myelin sheath up to a wavelength from its boundaries) by the input power in the guided modes at those particular wavelengths (before the paranodal region) to obtain the normalized transmission.

Bends

Bends are generated by extruding a circular cross-section along a sinusoidal path. The cross-section is not exactly a circle, but a 26 sided polygon with the vertices lying on the corresponding radius (for the axon and the myelin sheath). All these vertices then follow the sinusoidal path to construct a bend. The path is discrete too, with a step size equal to 0.5 μm. The number of vertices, and the step-size is optimized taking into account the accuracy and the speed of the simulation. With a straight path generated this way (discrete step size along the length, and a 26 sided polygon resembling a circle as the cross-section), and an eigenmode of the straight structure as the input, we ensure that we get close to 100% transmission. The percentage transmission is calculated by integrating the real part of the Poynting vector (with the fields corresponding to the output light) across the myelin sheath up to a wavelength away from the myelin sheath boundaries and dividing it by the source power. This is to include the evanescent fields and a few guided modes which are very close to the myelin sheath. Note that unlike the paranodal regions, we do not expand the output light in the basis of the guided modes at the end of the axon segment since the structure is continuously varying, and so is the basis of the guided modes. Some fraction of light in the non-guided modes at a particular cross-section might be included in the the basis of guided modes at an adjoining cross-section and vice-versa. Therefore it is more appropriate to observe the total power transmission instead of the modal transmission in such cases. We continue to be cautious, and ignore all the power inside the axon (a wavelength away from the myelin sheath boundary).

To account for the difference in the mode-profiles at different wavelengths, we send in the eigenmode corresponding to the central permissible frequency in a uniform straight axon, and observe the transmission in the myelin sheath up to a wavelength. We observe that for the wavelength corresponding to the central permissible frequency and lower, the transmission is close to 100%. But for longer wavelengths, the transmission can be substantially lower (e.g. for the thickest axon in our examples, the transmission is 96.81% for the longest wavelength). If the right mode (corresponding to the longer wavelength) had been incident, we would have obtained 100% transmission. To compensate for this insertion loss, we divide the transmission of the longer wavelengths by the transmission we obtain (for the same long wavelengths) when we send in a mode corresponding to the central frequency in a straight waveguide.

On rare occasions, for very small bends, this normalization procedure can yield slightly greater than 100% transmission (the maximum observed overshoot was ~0.18%) due to the finite resolution of the simulations, including the coarseness in the construction of the waveguide, and the import and export of field profiles across different programs (with different mesh sizes). In these cases the transmission is taken to be 100%. We also adopt this approach for other inhomogeneities which face this overshoot problem. We ran separate simulations with the exact input mode profiles for the particular wavelengths for a number of cases exhibiting the overshoot issue to verify that the transmission is indeed very close to unity in these cases.

We verify that the change of curvature seems to be the most important loss factor in the case of bends by running a few simulations for longer axonal segments (150 μm). For the same Δκ, the transmission for the longer segments was comparable to the transmission for the shorter ones.

Variable cross sectional area

The cross-sectional area of the myelin sheath is varied randomly according to an approximate normal distribution. We first generate 11 random points along a 50 segment with the desired mean and the standard deviation (s.d.) in Mathematica (assuming a correlation length ~5 μm–10 μm). As an example, for an axon with r = 3 μm, d = 2 μm, and s.d. 10% of d, the mean and s.d. of the points generated are 5 μm and 0.2 μm respectively. We fit these points with a polynomial of degree 7 (optimized over several trials). A polynomial of order 10 fits all the points exactly, but the local extrema of the function usually extend outside the span of the points it connects, and thus it has a greater randomness than that of the generated points. To ensure that the points in the fitted function indeed follow this distribution (with the expected mean and s.d.), we calculate the mean and the s.d. of this function by evaluating it at 200 points in the 50 segment. This process is repeated many times to get an appropriate function with the s.d. within 2.5% of the desired value. The Gaussian nature of the randomness of the function is ensured manually (by observing that ~95% of the points lie inside 2 s.d.). The normalization to account for the change of mode profiles with wavelength is exactly the same as that for the bends.

Non-circular cross-section

The non-circular cross-section (in the X-Y plane) is generated, similar in spirit to the varying cross-sectional area. Here, the random points, corresponding to the vertices of the cross-section of the axon, are generated separately for the 2 halves (one in the positive Y plane and one in the negative Y plane). In the positive Y plane, 10 points are generated at equal intervals from the polar angle 0 to n such that the mean separation of these points from the center is kept constant (equal to the mean radius). Now, a polynomial of order 7 is fitted to these points. The 2 points where this function crosses the X-axis (corresponding to the polar angles 0 and n) are noted. In the negative Y plane, 8 random points are generated at different polar angles. The other 2 points are those where the former function crossed the X axis. Now, a second fitting function (polynomial of degree 7) is generated with the weights of the couple of points lying on the X axis kept high to ensure that the function passes through these points. This is required because we want a continuity in the cross-sectional boundary for both the halves. This is the procedure for the construction of the cross-section of the axon. For the myelin sheath, we need to generate a parallel curve ensheathing the axon at a fixed perpendicular distance from the boundary of the axon. But a unique parallel curve for the g-ratio = 0.6 exists only when the s.d. of the boundary of the axon is small. For larger s.d. the segments in the generated parallel parametric curve start intersecting. Only an approximate parallel curve can be drawn in this case. We do that manually by selecting ~50 points separated from the axon’s boundary at the required perpendicular distance. Thus, the myelin sheath boundary is actually a ~50 sided polygon.

We choose a relatively long axonal segment (100 μm) and verify that almost all the non-guided modes of the waveguide are lost during propagation (by noting the transmission across many different cross-sections along the length and seeing that they converge). We integrate the real part of the Poynting vector just across the myelin sheath (not up to a wavelength) for each wavelength. We divide this transmission by the transmission just in the myelin sheath for the corresponding wavelength in a straight cylindrical waveguide of the same length, when the cylindrically symmetric mode (eigen mode for a circular cross-section) corresponding to the central permissible frequency is incident. This gives us an approximate normalized transmission for each wavelength. Following the procedure adopted while dealing with the inhomogeneities discussed earlier (e.g. bends, and varying cross-sectional area), we could have constructed separate inner and outer parallel curves at a wavelength separation from the myelin sheath and integrated the real part of the Poynting vector across that region. However, such unique parallel curves do not exist for long wavelengths and large inhomogeneities, and drawing approximate curves manually would also yield only approximately correct transmission values. We have verified that the results obtained using both procedures almost match with each other for a number of trial cases (within ~2% of each other). Since the transmission under such an inhomogeneity (the cross-sectional shape remains the same) almost does not drop with further increase in axonal length, slight inaccuracies in the transmission values do not matter.

Procedure for estimating the attainable transmission

We considered several optical imperfections to estimate the attainable transmission over the total length of an axon. We exponentiate the transmission fraction (obtained in our simulations for short axonal segments) for the variable cross-sectional area, the nodal and the paranodal regions, and the cross-talk between axons the required number of times. We do not exponentiate the transmission fraction for bends and non-circular cross-sections. For bends, as discussed earlier, we believe that the transmission depends primarily on the change of curvature (irrespective of the total length). For non-circular cross-sections, all the loss can considered as coupling loss (propagation loss is negligible). We then multiply all these transmission fractions to obtain the net transmission over the total length of an axon.

Additional Information

How to cite this article: Kumar, S. et al. Possible existence of optical communication channels in the brain. Sci. Rep. 6, 36508; doi: 10.1038/srep36508 (2016).

Publisher’s note: Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.