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.¹ In industrialized countries, PD affects 0.3% of all people and 1% of people over age 60.⁶ 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).⁴ 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.⁵ 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.

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

τ₁y’₁ = −y₁ + f₁(β₁ + h)

τ₂y’₂ = −y₂ + f₂(w₁₂y₁ + w₃₂y₃ + w₄₂y₄ − w₅₂y₅ + γ + b₂)

τ₃y’₃ = −y₃ + f₃(w₂₃y₂ + w₄₃y₄ + b₃)

τ₄y’₄ = −y₄ + f₄(w₃₄y₃ + b₄)

τ₅y’₅ = −y₅ + f₅(w₄₅y₄ + b₅)

γ = −w₆₂(−w₁₆y₁ − w₅₆y₅ + w₄₆y₄ + b₆)

Note that w₂₃ represents the difference between the excitatory and inhibitory inputs from TC to CT5. Note also that w_jk > 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.⁵ 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:

We can break down this system into 3⁵ = 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).

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:

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, b₆, so we took two estimates: b₆ = 6 for the low estimate, and b₆ = 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 b₆. The outputs for the low b₆ are shown below:

Healthy Solution:

Below is shown the firing rate outputs using these weights.

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

Parkinsonian Solution:

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 w₂₃ and w₃₂, which represent the excitatory and inhibitory connections between TC and CT5. Forcing all correlating weights to be equal in healthy and parkinsonian solutions except w₂₃ and w₃₂, we found a healthy solution that could be forced into a parkinsonian state by only altering w₂₃ and w₃₂. 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 w₂₃ and w₃₂. The region plot for the low b₆ is shown below.

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.

2. 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.

3. 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.

References

1. Sveinbjornsdottir, S. (2016).The clinical symptoms of Parkinson’s disease. Journal of Neurochemistry, 139(1), 318-324. https://doi.org/10.1111/jnc.13691.
2. DeLong, M. R., & Wichmann, T. (2007). Circuits and circuit disorders of the basal ganglia. Archives of Neurology, 64(1), 20–24. https://doi.org/10.1001/archneur.64.1.20.
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. https://doi.org/10.1146/annurev.ne.09.030186.002041
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. https://doi.org/10.3389/fnana.2015.00005.
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). https://doi.org/10.1142/S0129065718500211.
6. de Lau, L. M. L., & Breteler, M. M. B. (2006). Epidemiology of Parkinson’s disease. The Lancet Neurology, 5(6), 525-535. https://doi.org/10.1016/S1474-4422(06)70471-9.
7. Rall, W. (1955). Experimental monosynaptic input-output relations in the mammalian spinal cord. Journal of Cellular and Comparative Physiology, 46(3), 413–437. https://doi.org/10.1002/jcp.1030460303
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. https://doi.org/10.1016/j.neuroscience.2011.06.049
9. Nambu, A., & Llinaś, R. (1994). Electrophysiology of globus pallidus neurons in vitro. Journal of Neurophysiology, 72(3), 1127–1139. https://doi.org/10.1152/jn.1994.72.3.1127
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. https://doi.org/10.1523/JNEUROSCI.1719-05.2005
11. Hikosaka, O. (2007). GABAergic output of the basal ganglia. Progress in Brain Research, 160, 209–226. https://doi.org/10.1016/S0079-6123(06)60012-5
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. https://doi.org/10.1007/s002210050696
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. https://doi.org/10.1152/jn.1994.72.2.507
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. https://doi.org/10.1523/JNEUROSCI.23-05-01916.2003
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. https://doi.org/10.1016/s0306-4522(96)00555-6
16. Nambu, A. (2007). Globus pallidus internal segment. Progress in Brain Research, 160, 135–150. https://doi.org/10.1016/S0079-6123(06)60008-3
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. https://doi.org/10.1016/j.jtbi.2008.12.018
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. https://doi.org/10.1162/jocn.2010.21534
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. https://doi.org/10.1152/jn.00576.2010
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. https://doi.org/10.1155/2011/413543