Dynamical principles underlying song degradation in birdsong neural circuit

Our previous work has shown that AFP could adjust synaptic conductance of RA network with its firing pattern changed from bursting to tonic firing, eventually leading to song degradation. However, it is difficult to explore the dynamic of RA network due to its high dimensionality and complexity. To address this issue, we use the mean field theory proposed by Wilten Nicola to reduce RA network to a system of switching ordinary differential equations. Numerical simulations prove that the mean field model could qualitatively and quantitatively describe the behavior of nucleus RA. As a result, we perform the non-smooth bifurcation analysis over the mean field model to understand the bifurcations that RA network displays. Then we determine that there is a subcritical-Andronov-Hopf bifurcation at a particular applied current, which causes the emergence of the bistability with current increased. The bistability consists of a stable node and a stable limit cycle, corresponding to tonic firing and bursting of the full network, respectively. With these properties, we can explain the dynamical principles of AFP’s adjustment on RA during song degradation. Furthermore, our results would effectively promote the combination of dynamical analysis and neural network modeling based on the mean field theory.


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 HVC. c Our model shows two different firing patterns with different initial conditions within the same stimulus. d The effect of AFP on RA synapse conductance obtained by Henry. e Examples of syllables from Vikram Gadagkar. Reproduced with permission from references: [10], copyright 2016 AAAS. f Schematic of mean field theory Fig. 1d), we still do not 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 (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).
For this paper, we explain the dynamical principles of AFP's adjustment using the mean field theory and non-smooth bifurcation analysis with our article arranged as follows. Section 1 introduces the birdsong neural circuit and the phenomenon that we wish to explain with the mean field theory. Section 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. Most importantly, we explore the bifurcation that RA network displays and use it to explain the dynamical principles of AFP's adjustment. 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 NMDAR and AMPAR, to study the effect of T. Here, we use g R A /g R A = 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 . 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.

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-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 >: Equations (34-37) are our mean field model, which is used for analyzing RA network behavior.

Numerical simulation of the original network and the mean field model
To better understand the dynamics of RA, we numerically simulated our model with different applied currents. 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 , 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 <W(t)> The subcritical Hopf bifurcation that occurs in the mean field model and cross-sections at I Hop f = 1700 pA. a A subcritical Hopf bifurcation point appears at I Hop f ≈ 1641 pA, which causes stable nodes (red line and referring to tonic firing) and unstable limit cycles (grey cycles). As the current increases, the unstable limit cycle grows and eventually touches the stable nonsmooth 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 nonsmooth 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) assumptions during the derivation of the mean field model. When the current, 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. 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 Numerical simulation of the original network and the mean field model with the specific current. a and b are two different firing patterns with the same current and different initial conditions. c At I app =1700 pA, both bursting and tonic firing appears from different initial points. d AFP can adjust RA network with an initial condition in the yellow region and change normal syllables into distorted ones limit cycle and the node, which is hopeful to explain the network behavior and AFP's adjustment.

The non-smooth bifurcation analysis 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 Hop f , the trajectory would transform to a limit cycle from any initial point. But at I = I Hop f , a subcritical Hopf bifurcation point appears, and it would produce stable nodes and unstable limit cycles for I > I Hop f . 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 with a unique current 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),  >1645 pA). b The trajectory of the mean field model converges to the limit cycle (shown in red) when I app ≤1716 pA and to a stable node (shown in blue) when I app >1716 pA 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 principles behind AFP's adjustment during song degradation
Our prior study has shown that AFP can change normal syllables into distorted ones by adjusting RA's synapse structure. Even with the same current, 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 with some specific currents, 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 Hop f , a subcritical Hopf bifurcation point appears and then produces stable nodes and unstable limit cycles for I > I Hop f . Hence, we have the bistable structure of the mean field model when 1641 pA ≤ I app ≤ 1722 pA . The bistability consists of a stable node and a stable limit cycle, corresponding to tonic firing and bursting of the full network, respectively. During song degradation, AFP adjust synaptic conductance of RA network, which is equivalent to changing the initial value of the mean field system. Thus, during song degradation, the changed initial value would make the system translate to a stable node (corresponding to tonic firing and distorted songs of the complete network) rather than a stable limit cycle (corresponding to bursting and normal songs), eventually leading to song degradation. In conclusion, our work explains the dynamical principles of AFP's adjustment during song degradation based on the mean field theory and non-smooth bifurcation analysis and 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. Figure 5 shows the numerical simulation of the original network and the phase space trajectory of the mean field model with different currents. 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 [21] and some new functional neuron models which confirm the biophysical functions of biological neurons [22][23][24]. For example, Xiaohan Zhang [21] developed an RNNbased 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 [25]. 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 [26,27]. Weijie Ye [28] has performed a bifurcation analysis on this mean-field model and observed the synchronization states with different periods. Besides, La Camera et al. [29,30] 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.