Dynamical Bridge between Song Degradation and Neural Plasticity

High dimensionality and complexity are the main diﬃculties of the study over network dynamics. Recently, Wilten Nicola proposed the mean ﬁeld theory to research the bifurcations that the full networks display. Here, we use his approach on the birdsong neural network. Our previous work has shown that AFP could adjust the synapse conductance of nucleus RA and change RA’s ﬁring patterns, eventually leading to song degradation. To understand the dynamical principle behind this, we use a technique to reduce the RA network to a mean ﬁeld model, in the form of a system of switching ordinary diﬀerential equations. Numerical results have shown that the mean ﬁeld equations can qualitatively and quantitatively describe the behavior of nucleus RA. Based on the non-smooth bifurcation analysis of the mean ﬁeld model, we determine that there is a subcritical-Andronov-Hopf bifurcation at a particular stimulation, which can explain the role of AFP during song degradation. The results indicate that we can see AFP’s adjustment in RA synapse conductance as the adjustment of initial value in the bistable system. More precisely, during song degradation, the mean ﬁeld system would transform to a stable node (corresponding to distorted songs) rather than a stable limit cycle (corresponding to normal songs).


Introduction
Due to high similarity with human vocal neural circuitry [1], the birdsong neural circuit has been intensively studied [2]. As shown in Fig. 1a, the circuit involves song premotor pathway from the nucleus High Vocal Center (HVC) to the robust nucleus of the archistriatum (RA) and anterior forebrain pathway (AFP) [3]. AFP is related to the neuroplasticity of RA, especially in young and deaf birds [4]. For example, the neuroplastic changes under distorted auditory feedback (DAF) would result in birdsong degradation when a healthy adult bird is suddenly deafened [5]. Experimental proof has shown that AFP plays a decisive role here since song degradation would disappear with lesions of AFP [6].
The working principles of AFP during song learning and song degradation attracted researchers' attention. For example, Brainard [4] has discussed the potential instructive and permissive functions of AFP in vocal plasticity. Besides, Henry et al. [7,8] built a detailed biophysical model to simulate these processes. They found that AFP adjusts synapses in RA through the time difference, ∆T, between its input and HVC's input to RA (see Fig. 1b). Although they obtained the result about how ∆T affects synapse conductance of RA (see Fig. 1d), we still don't know how the neuroplastic changes induce some experimental phenomena, such as song degradation. Since the duration of some syllables is shorter than ∆T, AFP seems to only adjust the synapse conductance of nucleus RA and seldom affect RA during birdsongs [8]. To study AFP's adjustment, we have constructed a complete framework of the birdsong neural circuit before [9] and found that bursting in RA corresponds to normal syllables while tonic firing to distorted syllables. The RA network model in  this paper shows both these two patterns with different initial conditions and the same current stimulation (see Fig. 1e). Our results and other experiments [11] indicate that once an adult bird is deafened, ∆T will quickly affect the initial condition of RA and lead to song degradation. Here, we mainly focus on the nonlinear dynamic rationale behind this and try to build a dynamical bridge between the physiological performance and neuroplastic changes in RA based on the mean field theory introduced by Wilten Nicola [12].
Wilten Nicola has rigorously derived a mean field approximation to a large network of two-dimensional neurons based on population density methods and a quasi-steady state approximation. The mean field model is a system of switching ordinary differential equations. It can qualitatively and quantitatively describe the bifurcations that the full networks display. Wilten Nicola also considered the derivation of heterogeneous networks [13] or the network with white noise inputs [14]. Using the developed bifurcation theory of non-smooth systems [15], we can study the forces of inhibition, adaptation, and excitation interact in biological networks, such as the birdsong neural circuit being studied now (see Fig. 1f).
Our article is arranged as follows. Sect. 1 introduces the birdsong neural circuit and the phenomenon that we wish to explain with the mean field theory. Sect. 2 sets up the original network of nucleus RA and Sect. 3 derives the corresponding mean field model, using population density methods and a quasi-steady state approximation. In Sect. 4, we perform numerical simulation of these two models and establish the correspondence between them. Then we investigate the dynamic of the original network based on the non-smooth bifurcation analysis of the mean field model. Finally, we conclude this paper by discussing the differences between our two models and the comparison with other models in Sect. 5.

Method of Biological network
Before, we have constructed a complete framework of the birdsong neural circuit to study AFP's adjustment [9]. But here, we only model RA nucleus to investigate its dynamical mechanisms and consider ∆T between inputs from HVC and AFP. Based on the electrophysiological data [7,8], We set RA as a single all-to-all coupled Izhikevich network with its parameters defined in Table. 1. The dimensional form of RA network is given by for i=1, 2, . . . N. This model exhibits comparable physiological phenomena to experiments and other models [8,9]. During regular singing, RA neurons show bursting at 1255 pA, and once the bird deafened, RA neurons turn to show tonic firing, resulting in song degradation, as shown in Fig. 1c. For the further derivation of mean field systems, we rewrite this model in a dimensionless form which is given by The dimensionless parameters are also defined in Table. 1. As the dimensionless form is just the intermediate of original network and mean field model, we would mainly show the parameters in dimensional form below.
Besides, Henry established a simple biophysical model of synaptic plasticity, including mixed receptors of NM-DAR and AMPAR, to study the effect of ∆T. Here we use ∆g RA /g RA = H(∆T ) to simulate the result, which is given by 3 Method of the mean field model

Population density methods
For each neuron, we define x i = (v i , w i ) and it is a 2D point moving with time in the phase space X, which consists of vectors x = ( v, w). Then there will be N moving points in phase space. Now we consider a region Ω in it with piecewise smooth boundary and set the proportion of neurons in Ω as P (Ω, t). When N→∞ , we define the population density function, ρ(x, t), as follows where χ Ω is the conventional indicator function. The time evolution of ρ is a conservation equation, also called the continuity equation, given by Here J(x, t) is the flux, a vector of two components. Using Eqs. (12)(13), we can see J(x, t) = (J V (v, w, s, t), . Besides, the network averaged instantaneous firing rate, R i (t) , can also be described as Using Eqs. (10)(11)(12)(13)(14), we get a coupled ODE-PDE system to describe the dynamic of RA network, given by The system is numerically solvable, but not analytically solvable. Thus Eqs. (15)(16) need to be further simplified.

Simplification of the continuity equation
In this section, we use the moment closure assumption [16] to simplify the PDE (15). To begin, the population density function in its conditional form is and substituted into Eq. (15) to obtain Here we use the first order moment closure method to assume that w|v = w and integrate Eq. (18) with respect to w over the entire phase space W, getting Then we consider the mean adaptation variable, which is defined as Applying Eq. (11)(12)(13) and changing the order of integration as needed then gives: Then, using ρ(v, w, t) = 0 on ∂W and the boundary condition, we rewrite the equation (21) as Assuming that w w jump , we apply a Taylor expansion to yield Here we drop all the higher order terms and b in Eq. (23) since b 1. Hence, we only consider the O(1) problem: Finally, based on Eqs. (19)(24)(16), we obtain the following dynamical system:

The quasi-steady state approximation
In this model, the longest timescale is the adaptation variable time constant and the shortest time scale is that of membrane. Since the membrane time constant is 1 in the dimensionless system (5)(6)(7)(8), and a −1 = τ w 1. Then we introduce "slow time"t = τ −1 w t and obtain Since τ s = O(τ W ),we apply a quasi-steady state approximation to the PDE, which entails assuming that the density ρ(v, t) reaches its steady state density ρ(v), rapidly relative to s and w : Then we can solve Eq. (31) and get Noting that the denominator inside the integrand of Eq. (32) is not strictly positive on the phase space. A sign change of the flux corresponds to a saddle-node bifurcation for the denominator. Based on the analysis of Touboul [17], the bifurcation occurs when Hence, we give the expression of mean firing rate In conclusion, based on the the quasi-steady state approximation and other assumptions as discussed above, we reduce the system (5-8) to the following system of switching ordinary differential equations for s, < w >: The Eqs. (34-37) are our mean field model, which is used for analyzing RA network behavior.

Model simulations
Here we show how our model works under different stimulations. As shown in Fig. 2a, a random neuron displays bursting with its parameters shown in blue, such as voltage, V (t), adaptation current, W (t), and synaptic conductance, g(t). And the variation of average synaptic conductance shown in orange indicated that almost all neurons of the original network display bursting. At the same time, parameters of the mean field model, g syn s(t) and < w(t) >, regularly move with time, suggesting there might be a limit cycle (and it could be stable, as shown in Fig. 2c) in phase space. Then we can see the limit cycle of the mean field model as bursting of the original network. However, the duration of a limit cycle is shorter than that of bursting cycle due to the assumptions during the derivation of the mean field model. When the stimulation, I app , increased from 1200 pA to 2000 pA, neurons turned to display tonic firing, and the limit cycle seemed to become a stable node in the phase space of the mean field model, as shown in Fig. 2b and 2d. (b) (c)  Fig. 3 The subcritical Hopf bifurcation that occurs in the mean field model and cross-sections at I Hopf = 1700 pA. (a) A subcritical Hopf bifurcation point appears at I Hopf ≈ 1641 pA, which causes stable nodes (red line and referring to tonic firing) and unstable limit cycles (grey cycles). As stimulation increases, the unstable limit cycle grows and eventually touches the stable non-smooth limit cycle (blue surface), which corresponds to bursting. (b) The cross-section at I app = 1700 pA with a stable node (red point), an unstable limit cycle (grey cycle), and a stable non-smooth limit cycle (green cycle). (c) At I app = 1700 pA, mean field model with initial condition at the green region produces normal syllables and otherwise, it produces distorted syllables (corresponding to red region).
Based on the numerical simulation, we can consider the mean field model of switching ODE's as an adequate quantitative and qualitative descriptor of the behavior of the original network. More accurately, bursting in original network corresponds to a limit cycle in the mean field model and tonic firing corresponds to a node. The validation of the limit cycle and the node would be completed by MATCONT below. But what attracts us most is the dynamical change between the limit cycle and the node, which is hopeful to explain the network behavior and AFP's adjustment.

The dynamic of the mean field model
Based on the numerical continuation software MAT-CONT and direct simulations, we showed the dynamic of the mean field system in Fig. 3a, with I app as a bifurcation parameter. For I < I Hopf , the trajectory would transform to a limit cycle from any initial point. But at I = I Hopf , a subcritical Hopf bifurcation point appears, and it would produce stable nodes and unstable limit cycles for I > I Hopf . With I app increased, the unstable limit cycle grows and eventually collides with the stable non-smooth limit cycle.
These dynamics produce a bistable structure at a unique stimulation between 1641 pA and 1722 pA. For example, at I app = 1700 pA, Fig. 3b shows a stable node (red point), an unstable limit cycle (grey cycle), and a stable non-smooth limit cycle (green cycle) at the same time. In Fig. 3b, the stable node corresponds to tonic firing, and the stable limit cycle corresponds to bursting. Furthermore, initial points can be categorized into two types by the feature of syllables. As shown in Fig. 3c, RA neurons with initial condition in green region produce normal syllables and otherwise, they would produce distorted syllables.

The dynamical mechanisms behind AFP's adjustment
Our prior study has shown that AFP can change normal syllables into distorted ones by adjusting RA's synapse structure. Even under the same stimulation, the adjusted RA network shows tonic firing and produces distorted syllables. As shown in Fig. 4a, RA original network showed bursting with I app = 1640 pA, W (0) = 1255, and g(0) = 0. But when AFP increased synapse conductance from 0 to 40 nS, RA neurons showed tonic firing and induced distorted syllables as a result (see Fig. 4b). This phenomenon indicated a potential bistability that is hard to understand if we only consider the whole network. But now, we can explain the dynamical principle behind it with the mean field model. Based on the above numerical simulations, we found that the mean field model has the bistability at some specific stimulations, such as I app =1700 pA. For example, as shown in Fig. 4c, the difference in initial condition causes different dynamical results (limit cycle or node), which refer to bursting and tonic firing. We can establish a correspondence between the original network and the mean field model with their differences considered in the discussion below. Hence, the subcritical Hopf bifurcation point of the mean field can explain AFP's adjustment and the activity of the RA network. In fact, AFP's adjustment only needs to change the initial synapse conductance in RA to induce distorted syllables. As shown in Fig. 4d, RA network with an initial condition in the yellow region can be disordered by AFP to produce distorted songs. In contrast, the green region corresponds to RA networks which is hard to be disordered.

Conclusions
Our previous work has simulated the activity of the RA nucleus in a complete framework of the birdsong neural circuit, and we found that AFP could change the firing patterns of RA neurons from bursting to tonic firing during song degradation, resulting in distorted songs. Hence, we need a tool to understand the dynamic of RA nucleus under AFP's adjustment during song degradation and the tool we choose here is the mean field theory introduced by Wilten Nicola. The mean field model is a switching ODE's system derived from a large network of two-dimensional neurons based on population density methods and a quasi-steady state approximation. It can quantitatively and qualitatively describe the behavior of the original network. For example, bursting in the original network corresponds to a limit cycle in the mean field model, while tonic firing corresponds to a node. Then we can explain the dynamic of the original network by performing a non-smooth bifurcation analysis of the mean field model.
Based on the numerical continuation software MAT-CONT and direct simulations, it's amazing to found that at I = I Hopf , a subcritical Hopf bifurcation point appears and then produces stable nodes and unstable limit cycles for I > I Hopf . Hence, we have the bistable structure of the mean field model when 1641 pA ≤ I app ≤ 1722 pA . The bistability makes that the adjustment in the initial condition could change the behavior of the network. It can explain why AFP can degrade birdsongs by adjusting the synapse conductance of RA nucleus. In conclusion, our work builds a dynamical bridge between song degradation and neural plasticity. It supports the role of the mean field system as an adequate quantitative and qualitative descriptor of the behavior of the complete network.

The difference between the original network and the mean field model
The mean field model is a simplified form of the original RA network. It only considers the mean adaptation variable, w(t) , and the proportion of ion channels open in the membrane of RA neurons, s(t). Based on the numerical simulations, we found the mean field model is useful to qualitatively and quantitatively describe the dynamic of RA nucleus. However, there were some certain errors. For example, as shown in Fig. 2a, the duration of a limit cycle in the mean field model is shorter than that of bursting cycles in RA network. The reason for that is that the synaptic and adaptation is not large enough to satisfy the requirements of the quasi-steady state approximation.
Besides, the critical moment of transformation in the mean field model, from a limit cycle to a stable node, is different from that in the actual network, from bursting to tonic firing. Fig. 5 shows the numerical simulation of the original network and the phase space trajectory of the mean field model under different stimulations. For the original network, the moment of transformation is I app = 1645 pA , while for the mean field model, the moment of transformation is I app = 1716 pA, which is larger. Therefore, further improvements should focus on solving these two problems.

Comparison with other models
Recently, Alonso et al. proposed a model for canary song production, which considers the global activity of the units in each neural nucleus [18][19][20]. Alonso et al. and we all focus on the dynamic of the whole network, with differences laying in the model and research goals. They embedded neural additive models in an architecture compatible with the song system, hoping to create an artificial model to produce realistic songs. This paper established the mean field model to find the dynamic behind experimental phenomena and used them to inspired modeling and experiments. Combining the neural additive models with our mean field theory may be more useful to build a dynamical bridge between neural activity and experimental phenomena. Besides, the mean field theory can also help to analyze some complex decision-making tasks. For example, Xiaohan Zhang [21] developed an RNN-based Actor-Critic framework, which is trained through reinforcement learning (RL) to solve two tasks analogous to the monkeys' decision-making tasks. If we simplify the model to the mean field model, it will be easy to understand how RL trains the framework.
There are some other mean field models which we compare with our work here. Montbrió et al. proposed a new approach to reduce the large-scale neural network of quadratic integrate-and-fire neurons (QIF) and established an exact mean-field model [22]. The QIF network model could also be derived as a two-dimensional dynamical system for the firing rate and the mean membrane potential under the Lorentzian ansatz [23,24]. Weijie Ye [25] has performed a bifurcation analysis on this mean-field model and observed the synchronization states with different periods. Besides, La Camera et al. [26,27] derived a mean-field model for an uncoupled network of linear integrate and fire neurons with synaptically filtered noise considered. They used the same adaptation model and the same differential equation for the first moment of the adaptation current as us. These mean field models work on finding the dynamic of the full network from different angles, and it would be interesting to compare them in the same neural network.