SIAM News Blog

Nonsmooth Dynamical Systems in Neuroscience

By Wilten Nicola and Sue Ann Campbell

Large-scale models of the human brain, which help researchers understand humans’ rich plethora of potential behaviours, consist of millions of individual neurons coupled into large-scale networks. The Blue Brain Project, a European initiative that utilizes state-of-the-art computing resources, models large circuits of the brain in exquisite detail [6]. The Semantic Pointer Architecture Unified Network (SPAUN), a large functional network of 2.5 million model neurons, can even solve problems from the standardized IQ test [2].

While existing simulations replicate certain behaviours or experimental observations and thus provide an important step in understanding brain function, exhaustive parameter searches, additional simulations, and costly computing resources are necessary for real insight. A functional understanding of the behaviour of these large-scale models is needed to comprehend the human brain’s internal dynamics.

Fortunately, the macroscopic behaviour of large networks of identical (or even heterogeneous) subunits is often low-dimensional. We may derive a model, commonly called a mean-field or firing rate model, for macroscopic behaviour from the large-scale model. The derivation of mean-field equations frequently leads to models with discontinuities [3, 7], which arise because—at the most fundamental level—neurons have two states: quiescence (off) and firing (on). Analysis of the mean-field model can yield predictions and insights into the large-scale network that simulations alone would not provide. Prediction of the emergence of network-induced behaviour such as bursting, a special type of neuronal oscillation, is one such example. We may also derive mean-field models from first principles using knowledge of the underlying system or experimental data [10]. This approach helped to create a large-scale brain model that focusses on biologically-realistic connectivity [4]. In the mean-field approach, firing rate models represent the behaviour of a whole population of neurons, and can also be used as simplified models for individual neurons [3].

Firing rate models have led to applications both within and beyond neuroscience. For example, we can use firing rate equations to design networks of more complicated model neurons that display arbitrary dynamics [8]. This is the fundamental idea behind the brain simulator SPAUN; it presents a powerful potential application in the field of neuromorphic research, which aims to create engineered physical systems that process information in a manner similar to the human brain. Furthermore, networks of nonsmooth firing rate equations have surged in popularity with the growth of machine learning. In particular, deep learning uses networks of piecewise-smooth continuous nonlinearities to train large networks to solve what were once intractable machine learning problems. Examples include AlphaGo’s recent victory over world Go champion Lee Sedol [9].

The presence of nonsmooth nonlinearities in the models is the common theme of these applications. The analysis of nonsmooth dynamical systems has recently undergone rapid expansion due to the prevalence of phenomena in mechanical and electrical systems most easily explained by the assumption of switching behaviour in the underlying model [1]. The time is ripe to use these tools to understand the nonsmooth dynamical systems that arise in neuroscience and related applications.

Nonsmooth Dynamical Systems

An example of a nonsmooth dynamical system is the simple impacting system given by the following model:

\[\dot{y} = -2\eta y - x + F(t)\tag1 \]

\[\dot{x} = y\tag2 \]

\[x(t^-) = b, \:\: \Rightarrow x(t^+) = b, \:\: y(t^+) = -ry(t^-),\tag3 \]

where \(F(t) = sin(\omega t)\) is a simple periodic forcing function for an otherwise linear oscillator. Each time the oscillator impacts the surface \(x = b,\) the velocity changes discontinuously from \(y(t^-\) to \(-r\:y(t^-)\). If we were to remove the impacting boundary, the resulting system could only display oscillations. However, an impacting boundary radically alters the dynamical repertoire that this system can display, and even allows the system to exhibit chaotic behaviours for the right parameter regimes (see Figure 1). The equations of motion for the position \(x(t)\) and velocity  \(\dot{x}(t)\) of the cart form a nonsmooth system when accounting for interactions with the boundary.

Figure 1. Impact oscillator. 1a. The simple impact oscillator is derived from a cart attached to a spring with an external forcing term F(t). The impact oscillator can display behaviours that are not present for the system without an impacting boundary, such as chaotic dynamics. 1b. An example attractor. 1c.Trajectory on the attractor shown in 1b. Image credit: Wilten Nicola.

In the aforementioned example, the differential equation describing the system is smooth, but the impulses cause a discontinuous solution of the dynamical system. A dynamical system can also be nonsmooth if the right-hand side of the system itself is discontinuous. Such systems are defined via a set of differential equations with smooth right-hand sides, \(F_i,\) each valid on a different subset, \(S_i,\) of the phase space. We represent this mathematically as


\dot{x} = F_i(x, \mu), \: \:  x \in S_i, \: \:  i= 1,...,n, \: \:  x \in \mathbb{R}^n, \mu \in \mathbb{R}^p,\tag4 \]

where we assume that the boundary \(\Sigma_{ij},\) between regions \(S_i\) and \(S_j,\) commonly known as the switching manifold or discontinuity boundary [1], is a smooth \((n-1)\) dimensional manifold. Unlike more standard differential equations, nonsmooth differential equations can display different degrees of “nonsmoothness.” We classify the systems by the right-hand side’s smoothness, or equivalently, the solutions’ smoothness as they pass through the switching manifold. In Filippov systems, the solutions are continuous but have discontinuous time derivatives; in piecewise smooth-continuous (PWSC) systems, the discontinuity arises in the second order or higher time derivatives [1].

Nonsmooth systems described by (4) often display novel and exotic dynamical behaviours and bifurcations that the underlying smooth systems (which are separated by the switching manifold) cannot display themselves. While the current literature classifying and describing these nonsmooth dynamical behaviours is extensive [1], much work remains to generate a general bifurcation theory of nonsmooth dynamical systems that is as applicable and illuminating as the theory for smooth ones.

Nonsmooth Dynamics and Bifurcations in Neural Models

Depending on the underlying model, the dynamics of a single neuron can be quite complicated. Thus, the analysis of a network of neurons’ dynamical repertoire quickly becomes intractable. However, we can formally derive a low-dimensional system of differential equations describing the network behaviour in the large network or thermodynamic limit [3, 7, 10]. Due to their low dimensionality, these equations—called a mean-field model or firing rate model—are often much more tractable than the original network of neurons. Derivation of such a model from the original network is usually predictive of the behaviour of a sufficiently large network of model neurons. The bifurcation that causes the neuron to change from quiescence (a steady state equilibrium) to spiking (an oscillatory behaviour) plays a key role in these low-dimensional approximations, and can be represented by firing rate function \(f(I),\) which gives the frequency of the oscillation as a function of the current input into the neuron. When these curves are often nonsmooth (see Figures 2a and 2c), the resulting mean field or firing rate model is a nonsmooth dynamical system.

Figure 2. Nonsmooth neural model examples. 2a. The firing rate curve for type-I excitability often has a square-root type discontinuity. This leads to mean -field equations that are piecewise smooth-continuous (PWSC). 2b. The resulting mean-field systems display novel nonsmooth bifurcations, such as a saddle node on invariant circle (SNIC) boundary equilibrium bifurcation where a nonsmooth fold bifurcation collides with a nonsmooth limit cycle. 2c. A Heaviside function often approximates the firing rate curve for the Wilson-Cowan equations. This transforms the Wilson-Cowan equations into a Filippov system. 2d. The Filippov system features an intersection of two switching manifolds, which creates a new equilibrium in the plane. This equilibrium corresponds to the inhibitory stabilized network state in the related biological neural network. Image credit: Wilten Nicola.

Consider a generic network of type I model neurons, i.e., neurons that undergo a saddle-node on invariant circle (SNIC) bifurcation in their transition from quiescence to firing. Under appropriate conditions, we can formally show that the resulting mean-field model is approximately given by:

\[\tau_s \dot{s} = -s + \lambda_s F(s, w)\tag5 \]

\[\tau_w \dot{w} = -w + \lambda_w F(s, w)\tag6 \]

  F(s, w)=\begin{cases}
    \sqrt{H(s,w)}, & \text{$H(s,w) > 0$}\\
    0 & \text{$H(s,w) < 0,$}

where \(H(s, w) = \alpha_{00} + \alpha_{10}s + \alpha_{01}w + \alpha_{20}s^2,\)  with \(\alpha_{ij}\) determined by the parameters of the original model neurons [7]. The variables \(s\) and \(w\) correspond to network average values of currents in the neurons. In particular, \(s\) corresponds to excitatory synaptic coupling and \(w\) to inhibitory self-coupling, commonly referred to as spike frequency adaptation.

The system of equations governed by (5)-(6) is PWSC. A pair of local, nonsmooth, codimension 2 bifurcation points largely determines the entire bifurcation structure of this system. These points correspond to the collision of a saddle-node and Hopf bifurcation equilibrium point with switching manifold, \(H (s, w) = 0,\) respectively. As such, they act as organizing centers from which smooth saddle-node and subcritical Andronov-Hopf bifurcations and several nonsmooth, codimension 1 bifurcations emerge. This system also exhibits interesting nonsmooth limit cycle bifurcations, such as when a limit cycle becomes homoclinic with a boundary equilibrium bifurcation point [7] (see Figure 2b).

We can also derive a nonsmooth mean-field system using a Heaviside step function for the nonlinearities in the classic Wilson-Cowan model [10], leading to the following equations [5]:

\[u' = -u + f(\alpha_0 + u + \alpha_1 v)\tag8 \]

\[\tau v' = -v + f(\beta + u + \beta_1 v)\tag9 \]

    1 & \text{$x > 0$}\\
    0, & \text{x < 0.}
  \end{cases} \tag{10}

Unlike system (5)-(6), which only has a single switching manifold, equations (8)-(9) have a pair of switching manifolds (\(\alpha_0 + u + \alpha_1 v = 0\)  and \(\beta_0 + u + \beta_1 v = 0\)). Here we interpret the variables  and  as the firing rates of two groups of neurons in the network (excitatory and inhibitory neurons, respectively).

System (8)-(9) is a Filippov system. In one of the first analyses of a system with multiple switching manifolds, [5] demonstrates the existence of a nonsmooth analogue of the Hopf bifurcation occurring at the intersection points of these two switching manifolds for critical values of time constant \(\tau = \tau_{HB}.\) Furthermore, the authors illustrate the existence of nonsmooth analogues to the homoclinic and SNIC bifurcations.

Besides demonstrating interesting and novel nonsmooth bifurcations, the analysis of systems (5)-(6) and (8)-(9) yields substantial insights into the behaviour of the underlying neuron networks. We derived system (5)-(6) from a network of model neurons; simulations of the large network model [7] directly demonstrated the nonsmooth bifurcations. For example, the emergence of a nonsmooth limit cycle in the mean field model predicts when bursting will emerge in the network. System (8)-(9) demonstrates a link between the pseudo-focus occurring at the intersection of two switching manifolds and the inhibitory stabilized network (ISN) state [5] (see Figures 2c and 2d). In this state, feedback inhibition balances strong recurrent excitation to maintain stability. As we vary the parameters in the model equations, the ISN state loses stability via the nonsmooth Hopf analogue. Thus, nonsmooth analysis of these systems can yield substantial insights into the functioning of networks of neurons while simultaneously expanding the theory of nonsmooth dynamical systems.


[1] Bernardo, M., Budd, C., Champneys, A.R., & Kowalczyk, P. (2008). Piecewise-smooth Dynamical Systems: Theory and Applications (Vol. 163). London, UK: Springer.
[2] Eliasmith, C., Stewart, T.C., Choo, X., Bekolay, T., DeWolf, T., Tang, Y., & Rasmussen, D. (2012). A large-scale model of the functioning brain. Science, 338, 1202-1205.
[3] Ermentrout, G.B., & Terman, D.H. (2010). Mathematical Foundations of Neuroscience (Vol. 35). New York, NY: Springer.
[4] Falcon, M.I., Jirsa, V., & Solodkin, A. (2016). A new neuroinformatics approach to personalized medicine in neurology: The Virtual Brain. Current Opinion in Neurology, 29, 429-436.
[5] Harris, J., & Ermentrout, B. (2015). Bifurcations in the Wilson-Cowan equations with nonsmooth firing rate. SIAM Journal on Applied Dynamical Systems, 14, 43-72.
[6] Markram, H., Meier, K., Lippert, T., Grillner, S., Frackowiak, R., Dehaene, S.,… Grant, S. (2011). Introducing the human brain project. Procedia Computer Science, 7, 39-42.
[7] Nicola, W., & Campbell, S.A. (2016). Nonsmooth bifurcations of mean field systems of two-dimensional integrate and fire neurons. SIAM Journal on Applied Dynamical Systems, 15(1), 391-439.
[8] Nicola, W., Tripp, B., & Scott, M. (2016). Obtaining arbitrary prescribed mean field dynamics for recurrently coupled networks of type-I spiking neurons with analytically determined weights. Frontiers in Computational Neuroscience, 10, 15.
[9] Silver, D., Huang, A., Maddison, C.J., Guez, A., Sifre, L., Van den Driessche, G.,…Hassabis, D. (2016). Mastering the game of Go with deep neural networks and tree search. Nature, 529, 484-489.
[10] Wilson, H.R., & Cowan, J.D. (1972). Excitatory and inhibitory interactions in localized populations of model neurons. Biophysical Journal, 12(1), 1-24.

Wilten Nicola is a postdoctoral researcher in the Department of Bioengineering at Imperial College London. He completed his Ph.D. in applied mathematics at the University of Waterloo in 2015. Sue Ann Campbell is a professor of applied mathematics at the University of Waterloo. Her research interests include dynamical systems and computational neuroscience.

blog comments powered by Disqus