Biochemical basis of Quantum-like neuronal dynamics

The nervous system is a complex dynamical system that incorporates higher order biology (e.g., multicellular architecture) and lower-order biology (e.g., intra cellular pathway) that can be modeled via classical laws such as reaction-diffusion models. Simple reaction-diffusion models of neuronal tissue are shown to support bio-chemical wave effects that are analogous to quantum phenomena. These phenomena include quantum-like superpositions and classical entanglement which will not be affected by decoherence n the wet and warm brain environment. These classical phenomena could enable quantum-like complexity of brain functions. Conventional reactiondiffusion models of biological tissues challenge the current quantum brain hypothesis and suggest that the brain should perhaps be thought of as a classical quantum-like system. Statement of Significance This manuscript introduces the notion of nonseparability (classical entanglement) in the case of biochemical waves in arrays of coupled axons. We use a linear reaction-diffusion model with cross diffusion to address nonseparability between degrees of freedom (along and across the axon array). Perturbation theory applied to a nonlinear model with quadratic nonlinearity is used to illustrate nonseparability between modes along the axons. This paper suggests that the brain should perhaps be thought of as a classical quantum-like system.


Introduction
It has been argued that the complexity of brain functions cannot be explained by classical physical and/or chemical theories but requires quantum phenomena such as quantum superposition and entanglement to be understood [1] . While decoherence of quantum superposition in noisy and warm biological media serves as an argument against a quantum hypothesis for the brain functions, evidence of quantum mechanical effects in biological processes suggests that quantum phenomena could play a role in neurobiology [2] . For instance, it has been hypothesized that the nuclear spin residing in biological molecules may function as qubit, a quantum bit of information [3] , however, spin relaxation may be too fast to enable long lifetime quantum entanglement [4] . Outside the realm of neuroscience, the field of quantum biology has emerged from recent observations suggesting that a number of biological phenomena including enzyme catalysis, olfaction, photosynthesis may result from quantum mechanical effects such as coherence, tunneling and, entanglement [5] . In contrast to the quantum hypothesis of brain functions, there is an alternative proposal in which the complex dynamics of the nervous system leads to classical physical and chemical phenomena that are analogous to quantum mechanics and can therefore augment its capability for storing and processing information. Nonlocality and nonseparability are two distinctive aspects of entanglement. While nonlocality is a unique attribute of quantum mechanics, nonseparability is not. The notion of classical "entanglement " i.e., local nonseparable coherent superposition, or equivalently, the notion of classical nonseparability [ 6 , 7 ] has attracted interest and has been observed in the physical fields of optics [ 8 , 9 ] and more recently acoustics [ 10 , 11 ]. Since nonlocality of quantum phenomena (such as the Einstein's so-called spooky action at a distance ) is not a requirement for quantum information processing, classical nonseparability may offer the advantage of quantum superpositions in terms of information scalability without the fragility of quantum coherence. Classically entangled compositional variables in biological tissue, such as biochemical waves, each carrying information, would be interdependent. The classical inseparability of these variables, robust against decoherence, would enable quantum-like parallelism with potential relevance to understand the complexity of brain functions. The possibility of achieving quantum-like behavior in biological tissues does not require new theories but relies on reaction-diffusion models à la Turing based on classical laws [12] .

Reaction-diffusion model and propagating wave solutions
The purpose of the greatly simplified reaction-diffusion model introduced here is to demonstrate the emergence of quantum-like phenomena in a classical biological system. This model incorporates higher order biology (e.g., multicellular architecture) and lower-order biology (e.g., intra cellular pathway). The prototypical model is composed of three (infinitely long) parallel one-dimensional structures that can support propagative chemical waves [13] . We refer to these onedimensional structures as chemical waveguides. This architecture is reminiscent of parallel axons in the cerebral cortex [14] . It is also assumed that the chemical waves propagating along each chemical waveguide may result from an intracellular pathway involving a feedback loop between two chemical component, 1 and 2, subsequently called reactants. Inositol (1,4,5)-triphosphate (IP 3 )-mediated Calcium (Ca + ) waves in neurons are examples of waves that form from such a feedback process [15] . Indeed, the endoplasmic reticulum (ER) forms inside the axon a continuous network of interconnected tubules, sheets and cisternae which enables Ca 2 + -induced Ca 2 + release and re-uptake mediated by IP3, locally to and from the cytosol [16] . A continuous ER network is therefore able to support long-distance Ca 2 + signaling through the cytosol along the axon [17] .
The waveguides may also interact through ephaptic transmission between axons [18] or synaptic active or passive diffusion of reactants. Ephaptic coupling may result from exchange of ions between cells or may result from local electric or magnetic fields. Experimental evidence for axonal gap junctions suggests the possibility of ion exchange between neurons [19] . However, myelin sheaths surrounding the myelinated axons, the most common type of axons, in the central nervous system, are believed to inhibit direct neuron-to-neuron coupling. In contrast, unmyelinated axons are more common in the peripheral nervous system (e.g., in Remak bundles [20] ). However, in both cases of myelinated and unmyelinated axons, Schwann cells form close relationships with the axons. Myelinating Schwann cells surround a single axon whereas, unmyelinated Schwann cells bundle several axons together. The Schwann cells keep axons from touching each other. Schwann cells are glial cells subject to Ca 2 + signaling [21] . While the myelin sheath may still inhibit ion transport between Schwann cell and the axon, astrocytes (another type of glial cells) contact the axon membrane at the unmyelinated nodes of Ranvier. Contact between axon and glial cell at axon nodes offers a possible molecular mechanism by which neuron and glia can interact [22] . There is evidence that astrocytes influence the communication between neurons [23] . In unmyelinated nodes, gap junctions may connect glial cells and neurons [ 24 , 25 ]. Endogenous electric or magnetic fields can affect neuronal function via ephaptic coupling [26] . For instance, astroglial biomagnetic fields associated with Ca 2 + transients could be implicated in neuron-glial ephaptic crosstalk [27] .
Finally, while self-diffusion of the reactants can occur along and across the waveguides (e.g., via ephaptic coupling), the cross-diffusion of one reactant along or between waveguides may be driven by the gradient of the other reactant, that is, the two reactants diffusion may be mediated by cotransporters [28] . Another possibility for cross-diffusion across membranes is transport mediated by ubiquitous enzymatic transmembrane proteins. For instance, a hypothetical biomolecular mechanism for cross-diffusion across membranes was proposed [29] . A transmembrane enzymatic structure supporting two active and one regulatory site can transform a reactant 2 into a substance 3 at one of the active sites and vice versa at the other active site. The active sites sit on both sides of the membrane. The transformation is catalyzed by a reactant 1 at the regulatory site on one of the sides of the transmembrane structure. This mechanism assumes constant concentration of substance 3 and random distribution/orientation of the enzymic structures within membranes. Considering a chain of discrete compartments separated by membranes containing enzymic structures, the rate of change of reactant 2 in one compartment was shown to be proportional to the Laplacian of the of reactant 1; the signature of cross-diffusion. Some of the charac- teristics of this hypothetical mechanism have been observed in Receptor Tyrosine Kinases (RTK) [30] . RTK are constituted of a ligand-binding extra-cellular domain, a transmembrane domain and an enzymatic (kinase) intracellular domains. RTK can stimulate enzymatic activity leading to the cleavage of Phosphatidylinositol Biphosphate (PIP2) into IP3 and diacylglycerol (DAG). In turn IP3 production increases intracellular Ca 2 + concentration. Fig. 1 illustrates schematically the geometry of the model, irrespectively of the underlying biological mechanisms that may be involved in the self-diffusion, cross-diffusion, and reaction kinetics processes.
Initially, we assume small excursions in the concentration of reactants which enables us to consider a linear reaction model. That assumption will be lifted in a later section. The linear mathematical model takes the form of a set of coupled ordinary differential equations (ODEs) given by: The reactant compositional variables are ( ) 1 , 2 with = , , labeling the waveguides. The terms on the left of the equal signs represent the local accumulation of reactants. The diffusion constants d and d ′ quantify the effect of one reactant on its own transfer or on the transfer of the other reactant along the chemical waveguides, respectively. x is the positional variable along the waveguides. The transport coefficients L and L ′ describe transport between waveguides. The sign of these diffusion and transport coefficients determine the direction of transport with respect to concentration gradients. r ′ and r represent the self-regulation rate and cross regulation rate (feedback loop) of the intracellular pathway. The opposite signs in front of r ensures oscillatory behavior in concentration within a chemical waveguide. It is essential to note that equations (1a-f) do not introduce any new physics nor new chemistry.
Equations (1a-f) belong to a class of partial Differential Equations (PDE) known as cross-diffusion equations. Cross-diffusion systems are strongly coupled parabolic equations used to model processes arising in cell biology, thermodynamics and ecology [ 31 , 32 ]. Cross-diffusion systems obey parabolic equations of the form: where F i are nonlinear functions of the solution vector ⇀ = { 1 , 2 } . The nonlinear functions G i represent reaction terms. The cross-diffusion problem is well-posed if the Jacobian matrix of the vector field ⇀ is positive definite [ 33 , 34 ]. In the case of the Shigesada, Kawasaki, and Teramoto (SKT) model of interacting species [35] the cross-diffusion system possesses reaction terms. They showed that the cross-diffusion system can be approximated by a reaction-linear diffusion system and cross-diffusion instability is equivalent to Turing's instability.
In that context, we can reformulate Eq. (1a -f) as the system of equations: where the functions 1 1 , 2 1 , 1 2 , and 2 2 are here constants and related to the self and cross-diffusion coefficients for the direction along the waveguides d and d ′ and the transport coefficients L and L ′ for the direction perpendicular to the waveguides. Note that in Eq. (1a -f), the direction perpendicular to the waveguides is discrete. Furthermore, the gradient operator ∇ is continuous in the direction of the waveguides (direction x ) and discrete in the direction perpendicular to the planar array of waveguides (direction labelled by waveguide index = , , ). The finite number of waveguides imposes a condition of zero flux in the direction perpendicular to the waveguides at waveguides a and c . Equations (3a-b) are a generalized form of the SKT model. The functions g 1 and g 2 reduce to linear reaction rates in Eqs. (1a-f). The set of Eqs.
(1a-f) is therefore a linearized version of a SKT cross-diffusion system. The existence and nature of traveling wave solutions in excitable media with linear or nonlinear cross-diffusion has been reviewed. [36] .
With this in mind, we are seeking wave-like solutions of our linearized reaction cross-diffusion model of the form: ( ) where ( ) 1 , 2 represents the amplitude of the waves. "i " is the positive square root of -1. k is a wave number ( = 2 with being the wavelength) and is the angular frequency, one can compactify equations (1a-f) into a matrix form: is the matrix that couples the three waveguides. This is effectively the discrete Laplacian in the direction perpendicular to the wave guides including zero flux boundary conditions applied to the waveguides (a) and (c). The coupling matrix is a representation of the axonto-axon interactions in the multicellular architecture. In the present model, the three parallel waveguides are arranged in a single plane. Other configurations (e.g., the waveguides are arranged on the vertices of an equilateral triangle) will be associated with other forms of the coupling matrix. The dynamical matrix is given by . ⊗ is the usual tensor product operator. Equation (4) separates the reaction and diffusion dynamics along the waveguides from the coupling between them. It separates the dynamics of the system into a behavior along the x degree of freedom and a behavior corresponding to the degree of freedom associated with the labelling of individual waveguides. In that form, we can use a method developed by Othmer and Scriven [37] to seek a solution to Eq. (4) . We choose a solution for the amplitude with n being the eigenvalue associated with e n ) and where 2 × 1 = represents the amplitude of reactants 1 and 2 in any of the waveguides. Equation (4) reduces to The three eigenvectors e n corresponding to the eigenvalues 1 = 0 , 2 = −1 , and 3 = −3 , are: These eigenvectors represent the relative phase and magnitude of compositional waves supported by the different waveguides. They are signatures of the collective behavior of the interacting waveguides in our prototypical model of parallel axons. The waves supported by the three waveguides are effectively phase-locked. These eigenvectors are analogous to the Orbital Angular Momentum (OAM) degree of freedom in quantum mechanics.
Equation (5) reduces to solving the simpler single waveguide eigenvalue/vector problem We can now use Turing's linear stability analysis of our crossdiffusion systems to adjust the system to achieve plane wave traveling behavior of the solutions [38] . The parabolic equation (5) is the sum of several operators representing the self-diffusion (SD), cross-diffusion (CD) along and across the waveguides as well as the reaction kinetics (RK): = . The eigenvalues of the self-diffusion operator alone are real negative (since n < 0) and twice degenerate: The eigenvalues of the crossoperator are real positive or negative: = ±( − ′ 2 + ′ ) . The eigenvalues of the reaction kinetics are complex, = ′ ∓ . For a given wave number, self-diffusion alone leads to time-decaying concentration excursions and to a uniform stable state. Cross-diffusion alone leads to decaying solutions (stable uniform state) or growing (unstable state) for a given k . The reaction kinetics alone gives rise to damped or amplified temporal oscillations.
The eigen frequencies 0 n of the complete operator, sum of the SD, CD and RK operators, are found to be complex Complex eigen frequencies indicate that the compositional waves are generally attenuated or grow as they propagate. 0 n needs to be real (positive or negative) for the system to support traveling waves.
For this we are assuming that the system is in the limiting case where the rate r is sufficiently large such that 2 − ( − ′ 2 + ′ ) 2 is positive. This condition restricts propagation to waves with wave vectors that are sufficiently small. This corresponds to a strong feedback between reactants 1 and 2. Choosing the self-diffusion coefficients, d and L and the self-promoting reaction rate r ′ to be negligibly small leads to real frequencies corresponding to compositional waves propagating along the waveguides without attenuation. We see that the cross-diffusion introduces dispersion (dependency of the frequency on the wave number).
In that case, the complex eigenvectors are The reactants 1 and 2 oscillate with respect to each other with a phase of 2 (i.e., = 2 ).
In summary, choosing the second OAM as an example, the complete solution for the compositional wave supported by the coupled chemical waveguides is given by Solutions corresponding to the two other OAM can be written in a similar manner. This model can be easily generalized to a larger number, N , of coupled chemical waveguides (i.e., bundles of N > 3 axons) and various forms of the coupling matrix M N × N as well as larger number of reactants. Since equation (2) is a linear differential equation, any linear combination of solutions =1 6 × 1 , =2 6 × 1 and =3 6 × 1 is also solution of the equation. In the next section, we show that if the system is driven externally by an oscillating stimulus, we can obtain linear combinations of solutions which are interdependent in a way analogous to local quantum entanglement.

Externally driven reaction diffusion model
We solve the set of externally driven coupled differential equations Here d, r ′ and L ′ of equation (1a -f) are again neglected in order to focus on traveling compositional waves. F P with = 1 , 2 , 3 , 4 , 5 , 6 are the real amplitudes of the stimulus on each waveguide driving the two reactants. ( = 0 ) is the Dirac delta function that ensures that the stimulus is applied on each waveguide at the same location = 0 . is the angular frequency of the stimulus which may be different from the eigenvalues 0 n . From a biological point of view, the stimulus may be associated with excitation of apical dendrite cells from which emerge the axons forming the parallel fiber.
We now seek solutions of Eqs. (11a -f) in the form of linear combinations of propagative waves: These solutions are phase-locked with the stimulus i.e., the response wave of Eq. (12) is synchronized to the periodic frequency of the driver [39] .
We note in this equation that the wave number takes on different values corresponding to the different OAMs, namely k 1 , k 2 and k 3 . The six unknowns of the problem are the amplitudes, Z 1, 2 , Y 1, 2 and W 1, 2 . The system of equation (11a -f) is solved at = 0 . After extensive algebraic manipulations, we find ) and with 1 , 3 , 5 = 1 + 3 + 5 , 2 , 4 , 6 = 2 + 4 + 6 , 1 , 5 = 1 − 5 , 2 , 6 = 2 − 6 and 1 , 5 , 3 = 1 + 5 − 2 3 , 2 , 6 , 4 = 2 + 6 − 2 4 . The resonant amplitudes Z 1, 2 , Y 1, 2 and W 1, 2 are interdependent. Equation (12) , therefore, represents a coherent superposition. By tuning the driving amplitudes, all six resonant amplitudes change simultaneously. Similarly, by changing the frequency of the driver, all resonant amplitudes change as well. Importantly for the analogy with quantum mechanics, the resonant amplitudes are complex quantities. That is, they carry a phase. Equation (12) and ( 13a -c) establish an analogy with nonseparable quantum systems. This analogy can be best seen if we change to a notation conventional to quantum mechanics. Let us use the usual bra-ket notation of quantum mechanics. First, we denote the OAM eigenvectors of Eq. (6) by the kets | e 1 ⟩, | e 2 ⟩, | e 3 ⟩. These three vectors form an orthonormal basis for a three-dimensional Hilbert space, H OAM , i.e., they obey the conditions ⟨ | ⟩ = 1 if = and 0 otherwise. Here the bra ⟨e I | represents the complex conjugate (since the e I 's are real vectors, complex numbers and their conjugates are the same). Second, we also refer to 1 , 2 , 3 as the set of kets | k 1 ⟩, | k 2 ⟩, | k 3 ⟩. This set forms another orthonormal basis for a three-dimensional Hilbert space, H k , associated with the direction of propagation (i.e., x degree of freedom). The orthonormality condition implies the property ⟨ | ⟩ = ∫ dx − = 1 if = and 0 otherwise. With this notation, and considering separately the reactants 1 and 2, we can rewrite Eq. (12) in the form: Equations (14a , b ) are analogous to nonseparable superposition of states in the product Hilbert space H OAM ⊗H k . The basis set for states in that 9-dimensional product space is | e 1 It is impossible to write Eqs ( 12a,b ) as products of a linear combination of the kets: | e 1 ⟩, | e 2 ⟩, | e 3 ⟩ and a linear combination of the kets: | k 1 ⟩, | k 2 ⟩, | k 3 ⟩. For example, one cannot write C 1 as the product where the 's and 's are complex coefficients. This characteristic is the signature of nonseparability. Therefore, Equations (14a , b ) refer to nonseparable waves with the OAM degrees of freedom classically entangled with the directional degrees of freedom. The OAM forms a subsystem which can be described as the superposition of three mutually orthogonal states analogous to a qutrit [40] , which we call here a chem-trit. The direction of propagation forms another subsystem expressible as a superposition of three other orthonormal discrete states. It is also a chem-trit analogous to another qutrit. Therefore, equations (14a , b ) are describing the classical entanglement of these two chemtrits.
Exploring Eqs. (14) further, if we limit the stimulus of the system to amplitudes such that 1 , 5 = 1 − 5 = 0 and 2 , 6 = 2 − 6 = 0 , then 1 = 0 and 2 = 0 . Equations (14a , b ) reduce to In this case, both the OAM Hilbert space and the Hilbert space associated with the direction of propagation reduce to 2-dimensional spaces with bases | e 1 ⟩, | e 3 ⟩ and | k 1 ⟩, | k 3 ⟩, respectively. The OAM and direction of propagation form qubit analogues. We may call these subsystems chem-bits. C 1 and C 2 span a product space which is now a 2 × 2 = 4dimensional space with basis | e 1 ⟩| k 1 ⟩, | e 1 ⟩| k 3 ⟩, | e 3 ⟩| k 1 ⟩, | e 3 ⟩| k 3 ⟩. The states given by Eq. (15a ,b) are examples of two chem-bits in a nonseparable superposition analogous to the so-called Bell states in quantum mechanics. Other combinations of Bell states can be formed by restricting the stimulus amplitudes in various ways.
By varying the amplitude of the stimulus, the linear model system introduced above is able to explore a sizable portion of the product Hilbert space. For a given neuronal architecture involving N parallel axons and a set of R reactants, this product space scales linearly with the number of chemical waveguides as R × N . This model suggests that multiple bits of chemical information can be encoded in coherent superpositions and when the superposition is nonseparable its components can be processed efficiently in a simultaneous manner. This model system suggests the possibility of achieving quantum-like behavior in information processing with its high level of parallelism in classical physical/chemical systems such as neuronal tissues.

Nonlinear reaction diffusion model
The model system introduced in the preceding section depends on the assumption of intracellular pathways describable by linear functions of concentration. The nervous system, like many biological systems, is in fact a nonlinear dynamical system. When considering nonlinear reaction rate functions, it is possible to draw some observations on the emergence of other forms of nonseparable biological waves in the limit of small nonlinear effects using perturbation theory. For this, we augment the set of equation (11a -f) with a nonlinear term quadratic in the compositional variables.
All variables and parameters are the same as in Eqs (11a -f) but s which represents the rate of the quadratic reaction. We solve this set of differential equations using multiple time scale perturbation theory [41] in the limit = where ɛ is small. We rewrite the compositions as polynomials in ɛ , that is, considering as an example a generic composition, C : with the multiple time scales 0 = , 1 = , 2 = 2 . The time and space derivatives take the forms: Rewriting every composition and derivative of composition in Eqs (16a -f) according to equations (17) and ( 18a ,b), gives terms to zeroth order, first order and second order in ɛ . The zeroth order terms in ɛ can be grouped into a set of equations to zeroth order in perturbation: Equations (19a -f) are completely isomorphic to Eqs. (11a -f) which solutions are analogous to Eq (12) .
The terms linear in ɛ can be grouped to form the set of equations to 1 st order in perturbation: Equations (20a -f) form a system of driven differential equations. However, in contrast to the zeroth order equations which are externally driven, the first order equations are driven by the solutions to the zeroth order equations. The solutions of Eqs. (20a -f) are the sum of homogeneous solutions (solutions of the equation without the quadratic terms on the right-hand side of the equal sign) and particular solutions (i.e., with the right-hand side drivers). The time derivatives involving zeroth order compositions (e.g., ( ) 1 , 0 1 ) will lead to secular terms (i.e., diverging terms) in the homogeneous solutions. In order to eliminate secular solutions, we assume that the compositions are independent of the time scale 1 , that is the terms 1 are effectively zero.
Higher order equations are not necessary to support our argument in this section but can be constructed by grouping terms proportional to ɛ 2 .
We now recognize that the parallel fibers have a finite length on the order of a few millimeters. In that case the waves, solutions of the zeroth order equations, that can be supported by the driven chemical waveguides are not characterized by a continuum of wave number, k , but give rise to a set of waves with discrete values of the wave number, namely , ′ , ′′ , … This is a common property of waves in finite size systems. The relationship between the discrete wave numbers and the angular frequency for a given OAM state are still given by the real part of Eq. (8) . The solutions to the zeroth order equations are therefore a superposition of waves given by Eq. (12) over all possible wave numbers associated with a given OAM: Since Z 1, 2 , Y 1, 2 and W 1, 2 are resonant amplitudes (see Eq. (13a -c)), to simplify the problem, we consider that only the two wave numbers, k n and ′ , with frequencies, 0 n ( k n ) and 0 ( ′ ) closest to the driver's frequency, , will make significant contributions to C 6 × 1, 0 . Furthermore, to make our problem analytically tractable, we also assume that the external driver's amplitudes, F , are such that 1 , 2 = 0 . We therefore limit the model to excitation of only two OAMs. Eq. (21) reduces to: In that case, each wave characterized by a specific OAM is limited to two discrete levels.
Denoting ( ) = , ( ′ ) = ′ , ( ) = , ( ′ ) = ′ , we can use equation (22) to calculate the right-hand side terms of Eqs. (20a -f). As an illustrative example, we write one of them ( We can rewrite this zeroth order driving term as The first and second quadratic terms in Eq. (24) lead to selfinteraction between two states with the same OAM. The third term contains cross terms between states of different OAM, e 1 and e 3 . In contrast with the linear model, we treat these two OAMs as the signature of two subsystems, each subsystem supporting two discrete wavenumber levels. These are equivalent forms of chem-bits. The cross terms enable us to explore the tensor product Hilbert space of the Hilbert spaces of these two subsystems. The Hilbert space, H 1 , of the subsystem "1 " is two dimensional with basis: | 1 ⟩ = 1 and | ′ 1 ⟩ = ′ 1 . The two dimensional Hilbert space, H 3 , of the subsystem "3 " has the basis | 3 ⟩ = 3 and | ′ 3 ⟩ = ′ 3 . The tensor product space, 13 = 1 ⊗ 3 has dimension 2 2 with the basis 1 3 ) . With that notation, and for illustrative example, the cross term (CT) of Eq. (24) becomes: ( ( ) 1 , 0 ) 2 = 1 1 1 + 1 ′ 1 2 + ′ 1 1 3 + ′ 1 ′ 1 4 . The cross terms in Eq. (24) enable the exploration of the tensor product space of the bipartite system composed of the two subsystems "1 " and "3 ", i.e. two chem.bits. To that effect, we are seeking solutions to Eqs. (20a -f) which are only driven by the cross terms of Eq. (24) . A general form for these particular solutions in the chemical waveguide = , , is The question that arises is that of the separability of Eqs. (25a ,b). Let us consider a single chemical component in a single waveguide, (b), expressed in the bipartite system product basis, namely: Can we rewrite this expression as the product?
. This condition is not generally satisfied unless the waves for the two OAMs "1 " and "3 " are degenerate with same characteristic frequencies and wave numbers. Outside this very special case, the first order chemical waves given by Eqs. (25a ,b) are nonseparable in the tensor product space, H 13 . These nonseparable states are defined in H 13 which dimension 2 2 is exponentially complex. In the case of a neuronal tissue with N parallel axons, nonlinear reaction rates scaling as a power of N , and provided that each OAM supports only two k -levels, the multipartite nonlinear system would admit chemical waves to first order that span an exponentially complex Hilbert product space of dimension 2 N . The amount of information that could be encoded or processed by even a few tens of parallel fibers would then be massive. From a biological point of view, the spatially extensive apical dendritic cells in the cerebrum could be envisaged as a mean of sensing the wave number states, k , of nonseparable states. A wave number k relates to the wavelength of a chemical wave and its identification necessitates spatial resolution along the parallel fibers. The apical dendrites branch along the parallel fibers and form a large number of synapses [14] . The extensive synaptic network between apical dendrites and parallel fibers may be able to probe spatial modulations of reactants.

Conclusion
The biological wave models of a neurological tissue which have been presented in this paper (more specifically a highly reduced model of parallel axons) are based on conventional reaction diffusion theories. These models incorporate higher order biology (e.g., multicellular architecture through axon-to-axon coupling), lower-order biology (e.g., intra cellular pathway involving a number of chemical reactants) and external chemical stimulation. We have demonstrated theoretically that the reaction diffusion model with linear intracellular reaction rates offer the ability to prepare and tune nonseparable superposition of chemical waves that can propagate along the axons. These waves are products of a part associated with directional degrees of freedom along the chemical waveguides (axon) and an orbital angular momentum (OAM) part associated with the coupling between the multiple axons. Nonseparable quantum-like Bell states are constructed as a superposition of these chemical waves, which cannot be factored as a product of a single directional part and an OAM part. Remarkably, we also find that the amplitude coefficients of the nonseparable superposition of product waves are complex quantities, that is they include a phase. By tuning these complex amplitudes by varying stimuli distribution between the axons and/or frequency, the model is able to navigate a sizeable portion of the Bell state's Hilbert space. The dimension of this Hilbert space scales linearly with the number of chemical waveguides, N , times the number of chemical reactants (here 2).
We also showed theoretically, in the case of nonlinear (quadratic) intracellular reaction rates, using multiple time scale perturbation theory, the existence of another type of nonseparable superpositions of chemical waves that span a Hilbert space with exponential scaling. This space's dimension scales as the N th power of the number of waves that can be driven effectively by the external chemical stimulus. We show that this nonlinear reaction diffusion model is analogous to a two-partite (subsystems defined by two different OAMs), two-level quantum system (levels defined by the possible waves that can be stimulated). The amplitudes appearing in the nonseparable superposition are also complex quantities dependent on the frequency and amplitude of the external stimulus. By tuning these complex amplitudes, the stimulus can make the model navigate the exponentially scaling Hilbert space.
It must be admitted that the biological models presented here are oversimplified. The model multicellular architecture (parallel chemi-cal waveguides) and intracellular pathway (feedback between two reactants) may be a pale representation of the complexity of neuronal circuitry in the cerebral cortex. Nonetheless, these simple models demonstrate, within plausible conventional reaction diffusion models of biological systems, the remarkable ability of describing bio-chemical phenomena analogous to quantum mechanics. These phenomena include coherent superposition of states and nonseparability (i.e., local or classical entanglement). The nonseparability of bio-chemical waves is robust against decoherence as it is not a quantum phenomenon. However, it could enable quantum-like parallelism in bio-chemical wave communication and processing of information with potential relevance to understanding the complexity of brain functions. The possibility of achieving quantum-like behavior in conventional reaction-diffusion models of biological tissues challenges the current quantum brain hypothesis and opens up new avenues for interpretation of real neurological data and/or design of new experiments.

Declaration of Competing Interest
None.