Modeling Neuronal Systems as an Open Quantum System

: We propose a physical model for neurons to describe how neurons interact with one another through the surrounding materials of neuronal cell bodies. We model the neuronal cell surroundings, include the dendrites, the axons and the synapses, as well as the surrounding glial cells, as a continuous distribution of oscillating modes inspired from the electric circuital picture of neuronal action potential. By analyzing the dynamics of this neuronal model by using the master equation approach of open quantum systems, we investigated the collective behavior of neurons. After applying stimulations to the neuronal system, the neuron collective state is activated and shows the action potential behavior. We ﬁnd that this model can generate random neuron–neuron interactions and is appropriate for describing the process of information transmission in the neuronal system, which may pave a potential route toward understanding the dynamics of nervous system.


Introduction
The nervous system is a very complex part of living objects. It detects environmental changes that affect the body, and then it cooperates with the endocrine system in order to respond to such events. It coordinates its movements and sensory information through signal transmission with different parts of its body [1]. Therefore, understanding the structure and dynamics of the nervous system is the most challenging research area not only in biology but also in the development of modern science and technology, especially for the development of deep learning methods in a machine learning framework such as the computer vision task, speech recognition, multi-task prediction and the autonomous vehicle, etc. [2][3][4][5][6]. It may also help us solve the complex nervous system and assist in constructing proper implementations of the AI models associated with the understanding of the human brain [7]. Meanwhile, investigating the nervous system can not only enable us to eventually understand the mechanisms of human brain function but also reveal the origin of human thought and consciousness [8,9].
The key point of neural modeling is the connection between biological evidence and mathematical descriptions by using physical principles. Physically, we can depict all the phenomena into evolutions of physical states for simplification. Each neuron is commonly described by two states, the active (firing) and inactive (silence) states, and denoted physically as a spin or a quantum information qubit. The mathematical analogy of this concept was first proposed by McCulloch and Pitts in 1943 using a logic circuit for describing the neural networks [10]. About a decade later, Hodgkin and Huxley proposed the action potential model to explain neuronal state activation through the experimental evidence of signal transmission from the neuronal membrane [11]. They use electric circuits to mimic the lipid bilayer of membranes and the different kinds of ion channels through the membrane. This model can successfully fit the experimental data of the squid giant axon and other neural systems [11][12][13][14][15][16]. With the development of signal transmission and the binary neuron spin state, the concept of a neural network for brain modeling was also developed at almost the same time. Proposed in 1949, Hebb developed a theory to describe the formation of synapses [17]. He found that when two joining neuron cells fire simultaneously, the connection between them strengthens, which generates the new synapses. The concept of learning and memory explaining how the brain learns after firing several times through the connection between neuron cells was developed [17]. In 1982, Hopfield implemented Hebb's idea to simulate a neural network with a stochastic model, which is now called the Hopfield model. This model shows how the brain has the ability to learn and memorize [18]. In the Hopfield model, neurons are modeling as binary units or spins. The collective neurons form the state flow following the minimum "energy" principle to reach the fixed points, which are saved in advance. This has become the basic framework of the artificial neural network that AI development has relied on in the past four decades.
In the Hopfield model, neuron cells are randomly coupled to one another. It is known nowadays that the neural signals are transmitted between neurons mainly via neurotransmitters. It is not clear how the neural signals transmit through the randomly coupled neuron cells. In fact, the interactions between the neurons are realized through the neural surrounding materials, including the neurotransmitters through the synapses, the electrical signal from the dendrites and the axons and the noise from the glial cells surrounding the neurons [19]. Signal transmission via the intracellular wave propagation on the membrane has also been proposed from some evidence [20][21][22][23][24]. It is natural to ask how the neural surrounding affects the signal transmission between neurons, and it may be more fundamental to microscopically generate the random neuron couplings through these surrounding effects. As a result, the neural surroundings can be taken as an environmental effect with respect to the neuronal system as in the nonequilibrium dynamics description of open quantum systems from which the dynamics of macroscopic collective phenomena of the neuron system can be depicted.
In the literature, there are some approaches describing the process of neuron dynamics from the nonequilibrium dynamics of open quantum systems, where the collective phenomena of neuron systems are simulated as the spin collective phenomenon in open quantum systems determined by the master equation [25][26][27][28][29]. In the action potential picture, RC circuits are used to simulate the neural signals through the neuronal axon membrane and various ion channels. However, some evidence of inductive features in neural systems was observed [30,31]. Thus, the neural circuit for describing the signal transmission should be a neural RLC circuit. Physically, RLC circuits are equivalent to a collection of coupled harmonic oscillators with all kind of different frequencies [32], which provides us the idea to microscopically model the neural environment effect in the neuronal system as a continuous distribution of harmonic oscillators. On the other hand, the action potential is mainly due to the inward and outward ion flows of sodium and potassium ions through their own channels crossing the membrane, resulting in the voltage difference between the inside and outside of the membrane. The emergence of action potential transmission in the neurons is due to the stimulation from the environment by the neurons coupling with the surroundings. We exert our neuron system stimulations through the surrounding environment and demonstrate the collective neuron states with the action potential behavior. In this work, we will focus the investigation on the collective behavior of the neuronal system.
The paper is organized as follows. In Section 2, we introduce our neuronal model to describe the interaction of neurons in the nervous system through the surrounding neural materials. The surrounding environment effects are modeled as a continuous distribution of harmonic oscillators with different frequencies. Then in Section 3, we derive the master equation which governs the dynamics of neurons, where the random couplings between neurons naturally emerge. In Section 4, we analyze in detail the collective behavior of neurons in the neuronal system. The collective neuronal equations of motion are solved from the master equation obtained in Section 3. We apply external signals to stimulate the neuronal system, then the collective neuron states show the action potential properties. The thermal effects are also taken into account for mimicking the environment of the neuronal system. A conclusion is made in Section 5.

Modeling the Neuronal System as an Open Quantum System
Inspired by previous neuron modeling, such as the perceptrons [33], the linear associative net [34][35][36][37] and the content addressable memory [38][39][40], Hopfield modeled the neuron system based on neuron statistical behavior [18]. He started with the binary neural states. The idea of a binary neural state came from McCullough and Pitts's model [10], which characterizes the neural feature by a "all-or-none" logical property. In Hopfield's model, the binary neuron states can be realized through the spin configuration of N spins at time t: where α denotes different spin configurations and there are 2 N various configurations. The single spin state σ α i can either be 1 or −1 to represent the neuron activating or inactivating, respectively. The dynamical evolution of the neuron states is determined through the coupling matrix J N×N . The state evolution from time t to time t can be described with a matrix form of the transformation: which follows the equation where the sign function Sgn{} assigns the results 1 or −1 to the spin configuration elements with the time variable t , and the threshold of the action potential for ith neuron is represented by the bias b i . The dynamical evolution leads the many-neuron state to a local minimum in the configuration space. Thus, the Hopfield model can be equivalent to a disordered Ising model in statistics, where neurons are modeled as spins and the neuron dynamics are equivalently determined by an effective disordered Ising Hamiltonian: The first term in the above Hamiltonian describes the neuron coupling in the neuronal system, where σ z i is the z-component of the spin Pauli matrix for representing two states, silence (inactivate) and firing (activate), of the ith neuron. The second term is an effective magnetic field, in response to the threshold of the action potential for firing. Such a Hamiltonian mimics the signal transmitting between the neurons, in which all neurons are connected to one another with randomly distributed couplings J ij . By defining the coupling through the concept from Hebb's learning theory, this neural model has the ability to learn and memorize what it has saved. More specifically, the coupling strengths in Hopfield model are defined as J ij = ∑ ij ξ i ξ j , which comes from the random variable ξ i as a quenched, independent variable with equal probability in 1 and −1 [41].
Our motivation is to find the microscopic picture of the random neuron-neuron interaction from the interactions between the neurons and their surroundings (environment). Physically, neuron dynamics are governed by the interaction of neuronal cell body with their surroundings. The surrounding environment consists of all materials surrounding the neuronal cell bodies, including axon, dendrite, synapses and the surrounding glial cells. On the other hand, the neuronal system transmits the electrical signal that one can measure (see Figure 1a). Hodgkin and Huxley used the cable model to explain the electric voltage change of the neuronal membrane through RC circuits. However, the circuit analogy of the neuronal system should contain not only the resistance and capacitors from the fundamental action potential model but also the inductances as observed in some neuronal experiments [30,31]. Consequently, the neuronal system with neural signal transmission among neurons can be taken as more reasonable RLC circuits modified from the action potential model, which is shown in Figure 1b. In general, any circuit consisting of a complicated combination of RLC circuits corresponds to a collection of coupling harmonic oscillators [32,42]. Explicitly, a RLC circuit is equivalent to a damping LC harmonic oscillator: which can be obtained from a simple circuit equation L is the charge in the circuit, I(t) = dq(t)/dt is the corresponding circuit current, γ = R/L is the damping rate, ω = 1/ √ LC is the circuit frequency, and f (t) is an external force induced by the external voltage E (t) applied to the circuit. Quantum mechanically, a damping LC harmonic oscillator can be obtained microscopically from a principal LC harmonic oscillator linearly coupling to a continuously distribution of many LC harmonic oscillators in the surrounding environment with different frequencies: as is shown by Feynman and Vernon in their path-integral influence functional theory [32,43], where the first term is the principal LC oscillating Hamiltonian, the second term represents the environment consisting of a continuous distribution of the surrounding LC oscillators. The last term is the coupling of the principal oscillator with the environment oscillators, which results in the damping dynamics.
As a result, we model the neuron cells by a collection of spins, and these neuronal spins can interact with one another through the surrounding materials. The surrounding materials can be modeled by a continuous distribution of harmonic oscillators with all kind of different oscillating frequencies, characterized by Figure 1c, as inspired from the above RLC circuit picture. In quantum mechanics, the collection of the harmonic oscillators with all possible different frequencies can be expressed by the Hamiltonian in the particle number representation: where the continuous distribution of environmental oscillating modes has been taken ∑ k → (ω)dω and (ω) is the density of state of all the oscillating modes in the environment. The operators a ω and a † ω are the annihilation and creation operators of the oscillating mode with the corresponding frequency ω. The interaction between neuron spins and the surrounding materials is provided by the interaction Hamiltonian: where V ik is the coupling strength between the neuron spin mode σ z i and the oscillating mode ω k . The spin raising and lowering operators of each neuron spin are defined as σ ± i = 1 2 (σ x i + iσ y i ), and N is the total neuron number in the neuronal system. The neuron Hamiltonian of the following: corresponds to a set of neuron spins in an effective magnetic field g. With the above model description, we obtain our neuronal system described by the following Hamiltonian: which contains three parts: the Hamiltonian of the neurons, and the Hamiltonian of all oscillating modes from the surroundings, as well as the couplings between them. We will show that it is the coupling between the neurons and their surrounding oscillating modes that is causing the neurons' connection to one another and resulting in the random neuron-neuron coupling.

The Master Equation of the Neuron Dynamics
The neuron dynamics are described by the non-equilibrium evolution of all neuron states in the neuronal system. Neurons can interact with one another through their surrounding environment, as described by the Hamiltonian shown in Equation (9). Quantum mechanically, the total state evolution of neurons plus the environment is determined by the total density matrix ρ tot (t), which carries all the state information of the neurons and their surroundings. It is governed by the Liouville-von Neumann equation [44]: where H is the total Hamiltonian. We focus on the state evolution of the neurons affected by the infinite number of oscillating modes in the environment considered in Equation (9). Thus, we should focus on the time evolution of the reduced density matrix ρ S (t), which is determined by tracing over all the environmental states from the total density matrix: . The reduced density matrix ρ S (t) describes the nonequilibrium dynamical evolution of all neuron states in the neuronal system. The equation of motion for the reduced density matrix ρ S (t) that determines such evolution is called the master equation, which is derived below.
In the interacting picture of quantum mechanics, the total state of neurons plus their environment in the neuronal system is defined by ρ tot (t) = e iH 0 t ρ tot (t)e −iH 0 t , where H 0 = H S + H E . An expansion solution of Equation (10) in the interacting picture can be written as follows: where the interacting Hamiltonian is given by H I (t) = e iH 0 t H SE e −iH 0 t , and ρ tot (0) is the initial state. Suppose that the initial state of the system (neurons) and the environment (surroundings) is a decoupled state ρ tot (0) = ρ S (0) ⊗ ρ B (0), and the environment state is in a thermal equilibrium state, which is ρ B (0) = ρ B (0) = 1 Z B e −βH E , where β = 1 k B T is the inverse temperature of the environment and Z B = Tr[e −βH E ] is the environmental partition function. Meanwhile, we assume that all the neurons and their surroundings are weakly coupled to one another so that the environmental state almost remains unchanged, namely ρ tot (t ) ρ S (t ) ⊗ ρ B (0) in Equation (11), which is called as the Born approximation. Since the neurons are weakly coupled to their environment, we can further use the Markov approximation by replacing ρ S (t ) in the above Born approximation by ρ S (t). After making such Born and Markov approximations, we obtain the trace over all the environmental states. It can be shown that Tr E [H I (t), ρ tot (0)] = 0. Then, we obtain the master equation for the reduced density matrix ρ S (t) of all the neuron states.
This is known as the Born-Markov master equation in open quantum systems [45]. Now, we apply the above master equation formulation to the Hamiltonian described by Equation (9). For the sake of simplicity, we assume that the coupling strength is a real value, V i (ω) = V * i (ω). In the interacting picture, the interaction Hamiltonian is as follows.
After completing the trace over the environmental states in Equation (12) and changing the formulation back into the Schrödinger picture, we have the master equation only with neuron spin degrees of freedom: which describes the dynamical evolution of all the neuron states in the neuronal system. In Equation (14), the first term at the right hand side is a unitary transformation for the system dynamics according to the renormalized system Hamiltonian H S (t) = H S + δH(t), where δH(t) is induced by the couplings between the neurons and their surroundings:  (14): where κ ij and κ ij describes the effects of environmentally induced dissipation and fluctuation, and λ ij and λ ij are the environmentally induced random neuron-neuron couplings. The function J ij (ω) is the spectral density of the environment: and (ω) is the density of states of the environmental oscillating spectrum. The spectral density encompasses all the information about the structure of the materials surrounding neurons and the couplings with the neurons. The function n(ω, T) = Tr E [ ρ B (0)a † ω a ω ] is the particle distribution of the environmental oscillating modes andn(ω, T) = n(ω, T) + 1.
The neuronal system contains plenty of dynamical neurons and it is difficult to solve them even numerically. We can lower the calculating cost by summing up all the neuronal operators as a collective neural spin. The collective neural spinˆ S = (Ŝ x ,Ŝ y ,Ŝ z ) is defined by summing up all the neuron spins in each direction α: where α = x, y, x. In order to formulate the collective neuron behavior conveniently, we assume that the environment provides the same effect on all the neurons, namely, the coupling strength being independent relative to different neurons, V i (ω) = V(ω). The spectral density then becomes the following.
As a result, the master Equation (14) is simply reduced to the following form: where . The corresponding renormalization and dissipation/fluctuation coefficients in Equation (23) becomes the following.
This is the master equation for the collective neuron states in our physical modeling of neuronal system.

Equation of Motion for the Collective Neural States
The equations of motion for the collective neural states are obtained through the expectation values of the collective neural spin operators. By taking the expectation value of the collective neural spins with the reduced density matrix, as in the following: and applying the mean-field approximation, Ŝ αŜβ = Ŝ α Ŝ β , we obtain a close form of the equation of motion for the collective neural states from Equation (23): where the coefficients in the above equation of motion are provided by The coefficients K(t) and P(t) are related with the neuron-neuron interactions, and D(t) and F(t) are related to the dissipation and fluctuations. All of them are induced by the couplings between the neurons and with the environment in the neuronal system. In the following calculation, as an example, we obtain the most common spectral density as follows [43]: where η is a dimensionless coupling constant between the system and the environment, and ω c is a cut-off frequency. The value of s can be s = 1, < 1 and > 1, corresponding to the spectrum of the environment Ohmic, sub-Ohmic and super-Ohmic spectrum, respectively. Here, we consider the Ohmic spectrum, s = 1. With such detailed structure for the neuron environment, we can study the collective neuron behavior induced by the environment effects.

The Dynamics of Collective Neural States
With the equation of motion, Equations (29)-(31), obtained from the master equation, we want to explore how the collective neural states evolve in time under coupling with their surroundings. We start from the collective neural states (S x , S y , S x ) = (1, 1, −1) and study the time evolution of the collective state under different coupling constants. The results are presented in Figure 2, with the coupling strength being η = 0.1 and η = 0.02. As one can observe, the system always flows to the (0, 0, 0) state, which corresponds to the depolarized state. The differences are manifested in the trajectories of the collective neural states in the phase space. The larger the coupling constant, the sooner the collective neural state is depolarized. In reality, the neural signal transmissions through the external pulses stimulate the neurons. In the following, we want to explore how the collective neural states evolve in time when exerting an external pulse to the system by replacing the constant coupling strength with a rectangular pulse. In order to investigate the collective neural dynamics, we consider first the simple case of the nervous system at zero temperature and then move to the more realistic case of the neuron system at room temperature.

Collective Neural Dynamics at Zero Temperature
In order to consider an external pulse acting on the neuronal system, we modify the coupling as a time dependent parameter.
Meanwhile, at zero temperature, the dissipation and fluctuation coefficients in the master equation are reduced to the following: and κ(t) = λ(t) = 0 follows. Figure 3 shows the dynamics of the collective neural state after stimulation. We apply a simple square pulse with amplitude up to 0.8 for depolarizing the neural state (see the inset in Figure 3). The gray lines in Figure 3 represents time from 0 to 1.5 (ω 0 t), and the black-dashed line represents time from 1.5 to 10 (ω 0 t), where the duration time is scaled by the system frequency defined as ω 0 = 2000 Hz. Initially, the collective neural state is at polarized state (S x , S y , S z ) = (0, 0, −1). At the time from 0 to 1 (ω 0 t), no environmental noise exists, and the neuron system is in the rest state. At time from 1 to 1.5 (ω 0 t), we exert a pulse with amplitude 0.8, and the system reaches the depolarized state (S x , S y , S z ) ≈ (0, 0, 0). In the retrieving process, we extended the time duration to one and a half and made comparisons with the storing process and changed the coupling to 1.5 times less at time 1.5 to 2.25 (ω 0 t). The neural state then proceeds backward and experiences the repolarized process. Finally, during the time from 2.25 to 10 (ω 0 t), the system gradually returns to the initial rest state.  In Figure 4, we consider the stimulation when it doubles in time and becomes half in terms of quantity (refer to the inset of Figure 4). The result shows that the collective neural state can still repolarized to the original rest state. The similar action potential behavior was demonstrated from the collective neural state. This phenomenon shows that the stimulation of the neuronal system activates almost half of the neurons so that the collective neural state can reach the depolarized state "S Z 0". This is the "depolarization" mechanism of the neuron states via the environmental stimulation in our neuronal model.

Collective Neural Dynamics at Room Temperature
However, the zero temperature condition is an ideal case. The biosystem survives within room temperature. The increase in temperature will provide more noise from the environment. If we consider the environment at room temperature, all the dissipation and fluctuation coefficients of Equations (24)- (27) remain. At room temperature (T 300 K), the particle distribution in the environment can take on the classical Boltzmann distribution n(ω, T) = exp{−¯h ω k B T }. The state evolution under the room-temperature-environmental effect is shown in Figure 5. In this case, we find that it takes longer time for the state returning back to the initial rest state through the depolarization and repolarization processes due to the environment fluctuation on the collective neural state. This result also shows that no more than a half of the neurons are fired, and thus the maximum amplitude of the collective neural state is a little bit less than 0, but the temperature effect renders the firing states closer to the depolarized state. Furthermore, the condition of the same area also holds with the non-zero-temperature. The result is shown in Figure 6 for the same pulse profile in Figure 4.   Figure 4) in order to compare with the result in Figure 5. The room temperature is also the same T = 300 K.

Conclusions
In conclusion, we built a neuronal model as an open quantum system in order to understand the randomly coupled neuron-neuron interaction through their coupling with the neural environment. We used the master equation approach to studied the collective behavior of neurons incorporating the pulse stimulation in order to demonstrate the action potential. We explored the neuron dynamics at zero temperature and also at room temperature and found that, in both cases, the collective neural states evolved from polarized states (rest states) to the depolarized states and, finally, back to the initial polarized states under a simplified external pulse driving the neuronal system. Such results show that this simple neuronal model can not only catch the expected neuron dynamics but also provide an alternative mechanism in explaining how neurons couple or connect to one another through the complicated neuronal surrounding oscillating modes. As the result also shows, neuron-neuron connections through their surroundings are mainly determined by the spectral density of the neural surrounding environment, which characterizes the detailed energy spectrum structure of the neuronal environment. In this work, we only used the simple Ohmic spectral density as an example to simulate the collective neuron dynamics. The more realistic description of neuron dynamics relies on the spectral density that should be obtained from the spectral measurement of the neuron surroundings. Furthermore, a more complete description of the underlying neuron dynamics should be given by the neuron firing distribution in the neuronal system, which is depicted by the reduced density matrix of the master equation and can be obtained by solving the master equation directly. These remain open for the further investigations.