Model-based approaches to neuronal network firing and its subsequent validation with a previously recorded in-vivo dataset

The research featured in this blog post was performed by Carly Ferrell, Qile Jiang, and Olivia Leu, and the team was advised by Dr. Michael Caiola. This blog post was written by Carly Ferrell and published with minor edits. In addition to this post, the team has also given a midterm presentation, made a poster blitz video, and created a poster. They are currently working on a paper with the aim to publish it in an academic journal.

Mathematical Modeling of Healthy and Parkinsonian Firing Patterns in the Primate Thalamocortical Motor Circuit

Parkinson’s disease (PD) is a slowly progressing neuro-degenerative disease featuring impaired motor symptoms such as bradykinesia, muscular rigidity, and resting tremors.1 In industrialized countries, PD affects 0.3% of all people and 1% of people over age 60.6 The basal ganglia, motor thalamus, and motor cortex are three main components of the brain’s motor circuit and are responsible for movement planning and execution; movement disorders such as PD can develop when the typical activity of this circuit is disrupted.2,3 Specifically, PD is associated with the loss of dopaminergic neurons and altered neuronal oscillations in the beta-band (13-30 Hz).4 Other projects, such as the 2019 paper by M. Caiola and M. Holmes, have investigated the changes in the basal ganglia neuronal activity from a mathematical modeling perspective, but little research has been done on the parkinsonism-associated changes in the areas of the thalamus and cortex which are involved in the motor circuit.5 We employ a mathematical model to investigate network connection changes within the thalamocortical motor ciruit to better understand the transition from healthy to parkinsonian states in the brain.

Firing Rate Model

We choose to use a firing rate model to describe our system. This approach can successfully represent networks, since each unit in the model can represent a population of neurons receiving input (average firing rates) from other neuron populations.

A simplified circuit diagram of the thalamocortical motor circuit network is shown below, and provides the neuroscience basis for our model. The rounded squares each represent a population of neurons, which are connected by either excitatory (arrow-tipped lines) or inhibitory (circle-tipped lines) synaptic weights. The green circle represents the interneuron population of the thalamus.

Thalamocortical Loop Model

GPi (y1) globus pallidus internal
TC (y2) thalamocortical neurons
CT5 (y3) corticothalamic layer 5
CT6 (y4) corticothalamic layer 6
RTN (y5) thalamic reticular nucleus
IN (γ) thalamic interneuron population

Treating the interneuron population as a “relay,” γ, we can establish the following system of equations:

τ1y’1 = −y1 + f1(β1 + h)

τ2y’2 = −y2 + f2(w12y1 + w32y3 + w42y4w52y5 + γ + b2)

τ3y’3 = −y3 + f3(w23y2 + w43y4 + b3)

τ4y’4 = −y4 + f4(w34y3 + b4)

τ5y’5 = −y5 + f5(w45y4 + b5)

γ = −w62(−w16y1w56y5 + w46y4 + b6)

yi average neuronal population firing rate
wjk weight of the connection between populations j and k
h constant basal ganglia input
τi membrane time constant
fi activation function

Note that w23 represents the difference between the excitatory and inhibitory inputs from TC to CT5. Note also that wjk > 0 and τi > 0.

This can be represented with vectors and matrices as:

Ty’ = −y + F(x) ⟹ Ty’ = Ay + B

Activation Function Selection

Neurons traditionally respond to inputs sigmoidally.7,8,9 However, this model creates a nonlinear system of equations for which it is impossible to solve for eigenvalues analytically. In order to attain eigenvalues and be able to comment on the behavior of the model as a whole, we must establish a simpler activation function that still manages to approximate experimental neuron discharge behavior.5 A piecewise linear (PWL) activation function is ideal in our case, as it allows us to break down a complex system into linear pieces which can be solved and manipulated:

Piecewise Linear Activation Function

We can break down this system into 35 = 243 distinct linear regions in space, each with its own steady state (fixed point in space which the solution tends to as time increases). Out of these 243 regions, only the region in which each activation function is between 0 spikes/sec and its maximum firing rate contains a physiologically realistic steady state, further denoted as the middle region (outlined in green in the diagram below).

Accuracy of Piecewise Linear Activation Function to the Sigmoidal Activation Function

Data Matching

This semi-linear firing rate model has a number of constant values that we must locate in experimental data and incorporate, namely the baseline firing rates, maximum firing rates, and membrane time constants for each neuron population involved in our simplified motor circuit model. We were able to find values for these parameters through literature review, although some required that we make estimates informed by information from areas of the brain that behave similarly or data on these parameters from mice, rats, or cats. However, there does not seem to be data that documents the baseline firing rate for the thalamic interneuron population in the primate brain. Given our uncertainty about the true baseline firing rate value for the primate thalamic interneuron population, we decided to create two models, one with the low and one with the high baseline. The parameter values are shown in the table below:

Neuron Population bi Mi τi
GPi (y1) 55 Hz5,10,11,12,13 200 Hz14 8 ms11,15,16
TC (y2) 18.5 Hz17 300 Hz 25 ms
CT5 (y3) 7.25 Hz18 200 Hz 20 ms
CT6 (y4) 7.25 Hz18 200 Hz 15 ms
RTN (y5) 25 Hz17 500 Hz17 16.51 ms
IN (γ) Low: 6 Hz19
High: 22.7 Hz20

Stability and Steady State Conditions

No matter the disease state of our model, the neurons should not be at a state of maximal firing or absent firing for an extended period of time. Additionally, in Parkinsonian solutions, we should expect oscillations of firing rates. Thus the following must hold:

  1. Middle region contains its own steady state, and trajectories must not stabilize in another region.
  2. Healthy: Middle region is stable ⟶ trajectories are thus forced to stabilize in the middle region, making the system globally asymptotically stable.
    Parkinsonian: Middle region is unstable ⟶ trajectories are thus forced to oscillate around the middle region, forming a globally stable limit cycle.

To determine stability, the PWL activation function allows us to solve for the eigenvalues of each of the 243 regions explicitly. We found 3 possible cases:

  1. The region is stable regardless of weights.
  2. The region’s stability is conditional on weight values.
  3. The region (including the middle region) has eigenvalues that cannot be solved for analytically. Therefore, we used the Routh-Hurwitz Stability Criterion (RH) to derive 3 stability conditions.

Weight Search

The current literature does not specify the baseline firing rate for the interneuron population, b6, so we took two estimates: b6 = 6 for the low estimate, and b6 = 22.7 for the high estimate.

Comparing our data to the predicted values our model outputted, we were able to minimize the sum of squared error between the two and find a healthy solution for both the low and the high estimates of b6. The outputs for the low b6 are shown below:

Healthy Solution:

w12 = 1.520384442 w16 = 1.621278311 w23 = 0.4962387866 w32 = 1.117631687
w34 = 0.1540248925 w42 = 1.217895798 w43 = 0.0672671083 w45 = 1.542582263
w46 = 9.049109867 w52 = 4.5350845 w56 = 0.3330689302 w62 = 7.127373038

Below is shown the firing rate outputs using these weights.

Healthy Low FR

Here, all firing rates tend toward a specifc value as time increases, so they are stable solutions.

Parkinsonian Solution:

w12 = 1.520384442 w16 = 1.621278311 w23 = 0.8691494663 w32 = 0.5043792396
w34 = 0.1540248925 w42 = 1.217895798 w43 = 0.0672671083 w45 = 1.542582263
w46 = 9.049109867 w52 = 4.5350845 w56 = 0.3330689302 w62 = 7.127373038

Below is shown the firing rate outputs using these weights.


Here, several firing rates oscillate, indicating a limit cycle solution.

Weight Space

We were interested in the role of the thalamus in parkinsonian dysfunction, so we explored the relationship between w23 and w32, which represent the excitatory and inhibitory connections between TC and CT5. Forcing all correlating weights to be equal in healthy and parkinsonian solutions except w23 and w32, we found a healthy solution that could be forced into a parkinsonian state by only altering w23 and w32. In a parkinsonian solution, at least one of the Routh-Hurwitz stability conditions must be broken. We examined which condition or combination of conditions is broken when the system moves from a healthy to a parkinsonian state for different values of w23 and w32. The region plot for the low b6 is shown below.

Region Plot Low b6

Conclusions and Future Directions

  • Our model can represent the average firing rates of healthy and parkinsonian states in the thalamocortical motor circuit.
  • We established stability and steady state conditions for the system to be healthy or parkinsonian.
  • We have found multiple sets of weights that both satisfy the conditions and match the neuronal firing patterns in our in-vivo primate dataset.
  • We discovered that changing only the connection strength between TC and CT5 can force the system from a healthy to a parkinsonian state.
  • Next steps:
    • Examine the transition from both healthy to parkinsonian and parkinsonian to healthy.
    • Explore methods of biologically validating our model by pharmacologically manipulating the weights of motor circuit network connections.

More About the Team

  1. Carly Ferrell is a rising senior at Mississippi State University majoring in mathematics and minoring in statistics and music with a concentration in voice. She is interested in utilizing her skills in applied mathematics and statistcs to research music, specifically music theory and sight singing. Outside class, she enjoys reading, dancing, singing, and composing music.

Carly Picture

  1. Qile Jiang is a rising junior at Brown University majoring in Applied Mathematics. His primary research area is in applied dynamical systems, but he also has a keen interest in pure math topics such as algebra. Outside of school, he spends his time training boxing, painting, and going to operas and classical concerts.

Qile Picture

  1. Margaret Olivia Leu is a junior at Pomona College double majoring in mathematics and politics. She is interested in working on ways to use mathematics as a tool in the fields of politics and social justice work, and hopes to pursue a career that combines these two interests. Outside academics, she enjoys crocheting, cooking, and listening to music.

Olivia Picture


  1. Sveinbjornsdottir, S. (2016).The clinical symptoms of Parkinson’s disease. Journal of Neurochemistry, 139(1), 318-324.
  2. DeLong, M. R., & Wichmann, T. (2007). Circuits and circuit disorders of the basal ganglia. Archives of Neurology, 64(1), 20–24.
  3. Alexander, G. E., DeLong, M.R., & Strick, P.L. (1986). Parallel Organization of functionally segregated circuits linking basal ganglia and cortex. Annual Review of Neuroscience, 9(1), 357-381.
  4. Galvan, A., Devergnas, A., & Wichmann, T. (2015). Alterations in neuronal activity in basal ganglia-thalamocortical circuits in the parkinsonian state. Frontiers in Neuroanatomy, 9, 5.
  5. Caiola, M., & Holmes, M. H. (2019). Model and analysis for the onset of parkinsonian firing patterns in a simplified basal ganglia. International Journal of Neural Systems, 29(1).
  6. de Lau, L. M. L., & Breteler, M. M. B. (2006). Epidemiology of Parkinson’s disease. The Lancet Neurology, 5(6), 525-535.
  7. Rall, W. (1955). Experimental monosynaptic input-output relations in the mammalian spinal cord. Journal of Cellular and Comparative Physiology, 46(3), 413–437.
  8. Wilson, C. J., & Bevan, M. D. (2011). Intrinsic dynamics and synaptic inputs control the activity patterns of subthalamic nucleus neurons in health and in Parkinson’s disease. Neuroscience, 198, 54–68.
  9. Nambu, A., & Llinaś, R. (1994). Electrophysiology of globus pallidus neurons in vitro. Journal of Neurophysiology, 72(3), 1127–1139.
  10. Kita, H., Tachibana, Y., Nambu, A., & Chiken, S. (2005). Balance of Monosynaptic Excitatory and Disynaptic Inhibitory Responses of the Globus Pallidus Induced after Stimulation of the Subthalamic Nucleus in the Monkey. Journal of Neuroscience, 25(38), 8611–8619.
  11. Hikosaka, O. (2007). GABAergic output of the basal ganglia. Progress in Brain Research, 160, 209–226.
  12. Wichmann, T., Bergman, H., Starr, P. A., Subramanian, T., Watts, R. L., & DeLong, M. R. (1999). Comparison of MPTP-induced changes in spontaneous neuronal discharge in the internal pallidal segment and in the substantia nigra pars reticulata in primates. Experimental Brain Research, 125(4), 397–409.
  13. Bergman, H., Wichmann, T., Karmon, B., & DeLong, M. R. (1994). The primate subthalamic nucleus. II. Neuronal activity in the MPTP model of parkinsonism. Journal of Neurophysiology, 72(2), 507–520.
  14. Hashimoto, T., Elder, C. M., Okun, M. S., Patrick, S. K., & Vitek, J. L. (2003). Stimulation of the Subthalamic Nucleus Changes the Firing Pattern of Pallidal Neurons. The Journal of Neuroscience, 23(5), 1916–1923.
  15. Nakanishi, H., Tamura, A., Kawai, K., & Yamamoto, K. (1997). Electrophysiological studies of rat substantia nigra neurons in an in vitro slice preparation after middle cerebral artery occlusion. Neuroscience, 77(4), 1021–1028.
  16. Nambu, A. (2007). Globus pallidus internal segment. Progress in Brain Research, 160, 135–150.
  17. van Albada, S. J., & Robinson, P. A. (2009). Mean-field modeling of the basal ganglia-thalamocortical system. I Firing rates in healthy and parkinsonian states. Journal of Theoretical Biology, 257(4), 642–663.
  18. Opris, I., Hampson, R. E., Stanford, T. R., Gerhardt, G. A., & Deadwyler, S. A. (2011). Neural Activity in Frontal Cortical Cell Layers: Evidence for Columnar Sensorimotor Processing. Journal of Cognitive Neuroscience, 23(6), 1507–1521.
  19. Ison, M. J., Mormann, F., Cerf, M., Koch, C., Fried, I., & Quiroga, R. Q. (2011). Selectivity of pyramidal cells and interneurons in the human medial temporal lobe. Journal of Neurophysiology, 106(4), 1713–1721.
  20. Putrino, D. F., Chen, Z., Ghosh, S., & Brown, E. N. (2011). Motor Cortical Networks for Skilled Movements Have Dynamic Properties That Are Related to Accurate Reaching. Neural Plasticity, 2011, 1–15.
Mike Caiola
Mike Caiola
Affiliate Scientist, Mentor