Analysis and Control of the Dopamine Circadian Rhythms Model

Lakshmi N Sridhar; Lakshmi N Sridhar

ISSN: 2455-5282

Global Journal of Medical and Clinical Review Articles

Review Article Open Access Peer-Reviewed

Analysis and Control of the Dopamine Circadian Rhythms Model

Lakshmi N Sridhar*

Chemical Engineering Department, University of Puerto Rico, Mayaguez PR 00681-9046, Puerto Rico

Author and article information

*Corresponding author: Lakshmi N Sridhar, Chemical Engineering Department, University of Puerto Rico, Mayaguez PR 00681-9046, Puerto Rico, E-mail: [email protected]

doi : 10.17352/2455-5282.000195

Received: 23 January, 2025 |Accepted: 01 February, 2025 | Published: 03 February, 2025

Keywords: Circadian; Dopamine; Bifurcation; Optimal control

Cite this as

Sridhar LN. Analysis and Control of the Dopamine Circadian Rhythms Model. Glob J Medical Clin Case Rep. 2025:12(2):029-036. Available from: 10.17352/2455-5282.000195

Copyright License

© 2025 Sridhar LN. This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

Abstract

Background: The high nonlinearity of the dopamine circadian rhythms model is seen in the presence of limited cycles that disrupt the circadian rhythms. Limit Cycles originate from Hopf bifurcation points. Bifurcation analysis and Multiobjective nonlinear model predictive control are performed on the dopamine circadian rhythms model.

Methods: The MATLAB software MATCONT was used to perform the bifurcation analysis. The Multi-objective Nonlinear Model Predictive Control was performed using the optimization language PYOMO.

Results: The Bifurcation analysis reveals Hopf Bifurcation points that produce limit cycles. To eliminate the rhythm disturbing limit cycles the bifurcation parameter is multiplied by an activation factor involving the tanh function. The nonlinearity of the dopamine circadian rhythms model also causes spikes in the control profiles when multiobjective nonlinear model predictive control calculations are performed. The spikes are also eliminated when the control variable is multiplied by the same activation factor.

Conclusion: The dopamine circadian rhythms model is shown to have two Hopf bifurcations, which cause limit cycles that can disrupt the circadian rhythms. An activation factor involving the tanh function eliminates the limit cycle causing Hopf bifurcations. This activation factor also removes the spikes that occur in the control profile.

Main article text

Background and introduction

Graybiel and Kimura [1] researched the basal ganglia in adaptive motor control. Vitaterna, et al. [2] studied the differential regulation of mammalian period genes and circadian rhythmicity by cryptochromes 1 and 2. Van der Horst, et al. [3] demonstrated that mammalian cry1 and cry2 are essential for the maintenance of circadian rhythms. Strogatz [4] explored the onset of synchronization in coupled oscillators. Shearman, et al. [5] discussed how the molecular loops interact in the mammalian circadian clock. Abarca, et al. [6] showed that Cocaine sensitization and reward are under the influence of circadian genes and rhythm. Preitner, et al. [7] demonstrated that the orphan nuclear receptor rev-erbα controls circadian transcription within the positive limb of the mammalian circadian oscillator.

Lotharius and Brundin [8] investigated the pathogenesis of Parkinson’s disease with regard to dopamine vesicles and alpha-synuclein. Oleksiak, et al. [9] studied the variation in gene expression within and among natural populations. Boeuf, et al. [10] investigated the individual variation of adipose gene expression and identified covariate genes by cDNA microarrays. Forger and Peskin [11] provided a detailed predictive model of the mammalian circadian clock. Leloup and Goldbeter [12] came up with a detailed computational model for the mammalian circadian clock. Sato, et al. [13] demonstrated a functional genomics strategy that reveals rora as a component of the mammalian circadian clock. Castaneda, et al. [14] discussed the circadian rhythms of dopamine, glutamate, and GABA in the striatum and nucleus accumbens of the awake rat.

Kienast and Heinz [15] conducted research on Dopamine in the diseased brain. Sigal, et al. [16] discussed the Variability and memory of protein levels in human cells. Hong, et al. [17] showed a strategy for robust temperature compensation of circadian rhythms. McClung [18] researched the effect of circadian genes, rhythms, and the biology of mood disorders. Sleipness, et al. [19] discussed the diurnal differences in dopamine transporter and tyrosine hydroxylase levels in rat brain: Dependence on the suprachiasmatic nucleus. Hampp, et al. [20] investigated the regulation of monoamine oxidase by circadian clock components.

Liu, et al. [21] investigated the redundant function of rev-erbα and β and the non-essential role of bmal1 cycling in the transcriptional regulation of intracellular circadian rhythms. Hood, et al. [22] demonstrate that Endogenous dopamine regulates the rhythm of expression of the clock protein per2 in the rat dorsal striatum via daily activation of d2 dopamine receptors. Ripperger, et al. [23] investigated the daily rhythm of mice. Gravotta, et al. [24] showed that the depletion of dopamine using intracerebroventricular 6-hydroxydopamine injection disrupts normal circadian wheel-running patterns and period2 expression in the rat forebrain. Bugge, et al. [25] show that the Rev-erbα and rev-erbβ co-ordinately protect the circadian clock and normal metabolic function. Solt, et al. [26] demonstrated a method to regulate circadian behavior and metabolism by synthetic rev-erb agonists. While Ikeda, et al. [27] discovered the molecular mechanism that regulates the 24-hour rhythm of dopamine d3 receptor expression in mouse ventral striatum.

Karatsoreos [28] studied the links between circadian rhythms and psychiatric disease. Jager, et al. [29] investigated the behavioral changes and dopaminergic dysregulation in mice lacking the nuclear receptor. Chung, et al. [30] analyzed the impact of circadian nuclear receptor rev-erba on midbrain dopamine production and mood regulation. Ye, et al. [31] discuss bmal1 inhibition mediated by cryptochrome and period proteins in the mammalian circadian clock. Fifel, et al. [32] showed how the daily and circadian rhythms were altered following dopamine depletion in mptp-treated non-human primates. Colwell [33] studied the link between neural activity and molecular oscillations.

Bedrosian, et al. [34] investigated the timing of light exposure affects mood and brain circuits. Huang, et al. [35] showed that Circadian modulation of dopamine levels and dopaminergic neuron development contributes to attention deficiency and hyperactive behavior. Lee, et al. [36] identified a novel circadian clock modulator controlling bmal1 expression through a ror/rev-erb-response element-dependent mechanism. Takahashi [37] studied the transcriptional architecture of the mammalian circadian clock. Albrecht [38] studied the effect of molecular mechanisms in mood regulation involving the circadian clock. A comparison of macroscopic models for human circadian rhythms was conducted by Hannay, et al. [39]. Kim and Reed [40] developed a mathematical model involving circadian rhythms and dopamine.

This research aims to conduct a bifurcation analysis and perform multiobjective nonlinear model predictive control (MNLMPC) calculations on the dopamine circadian rhythms model of Kim and Reed [40]. The bifurcation analysis reveals the existence of unwanted Hopf bifurcations that cause limit cycles. An activation factor is used to eliminate the Hopf bifurcations. The MNLMPC calculations revealed spikes in the control profiles. These spikes are also eliminated using the same activation factor.

The paper is organized as follows. The dopamine circadian rhythms model of Kim and Reed [40] is first described. The bifurcation analysis and multiobjective nonlinear model predictive procedures are then described followed by the results and discussion. The conclusions are then presented.

Dopamine circadian rhythms model

In the Dopamine Circadian model by Kin and Reed [40] the variables are X = [P₁,P₂,P₃,P₄,C,PC,PC_N,PN,C_N,BC,REV,ROR,S,TH,MAO,DRD3]. P₁,P₂,P₃,P₄ represent the period proteins PER with 1 2 3 4 i phosphorylations in the cytosol. C represents the CRYPTOCHROME protein (CRY) in the cytosol. PC, PC,PCN,PN,CN represent the PER-CRY protein dimer in the cytosol, the nuclear PER-CRY protein dimer that inhibits BMAL1-CLOCK, the nuclear PER, and the nuclear CRY; that inhibits transcription driven by BMAL1-CLOCK. BC,REV,ROR,S,TH,MAO,DRD3 represent the BMAL1-CLOCK; a protein dimer that drives the transcription of core clock genes, the orphan nuclear receptor REV-ERB, the retinoic acid receptor-related orphan receptor, the BMAL1 circadian clock gene, Tyrosine hydroxylase which is a rate-limiting enzyme in dopamine synthesis, the Monoamine oxidase which catalyzes dopamine degradation and the Dopamine receptor.

Table 1 contains the description of the model variables. The model consists of 16 differential equations and has three linked parts. The first part is the core circadian clock that consists of a feedback loop containing BC, P₁,P₂,P₃,P₄ PC, PC,PC_N,PN.

The second part is the secondary feedback loop involving REV, ROR, S, and BC that influences the core clock while the third part is the dopamine (DA synthesis the elements of which are REV, ROR, and DRD3. All concentrations are in nanomoles (nM) and all rates are in nanomoles per hour (nM/hr). More details can be found in the work of Kim and Reed [40].

The model equations include algebraic expressions that are used in time-dependent differential equations. The function f is defined as

$f (x, y, z) = \frac{(x - y - z) - {[{(x - y - z)}^{2} + 4 z x]}^{1 / 2}}{2 x} (1)$

While the other functions are

$B C_{f r} = f (B C, P C_{N}, k_{d}) (2)$

$R_{s} (S, R E V) = \frac{ρ_{s}}{{(1 + k_{s} (1 - f (S, R E V, ε_{s})))}^{n_{s}}} (3)$

$A_{s} (S, R E V, R O R) = α_{s} f (S, R E V, ε_{s}) (\frac{R O R}{R O R + k_{s}}) (4)$

$R_{t h} (T H, R E V) = \frac{ρ_{t h}}{{(1 + k_{t h} (1 - f (S, R E V, ε_{t h})))}^{n_{t h}}} (5)$

$A_{t h} (T H, R E V, R O R) = α_{t h} f (T H, R E V, ε_{t h}) (\frac{R O R}{R O R + k_{t h}}) (6)$

$R_{d r} (D R D 3, R E V) = \frac{ρ_{d r}}{{(1 + k_{d r} (1 - f (D R D 3, R E V, ε_{d r})))}^{n_{d r}}} (7)$

$A_{d r} (D R D 3, R E V, R O R) = α_{d r} f (D R D 3, R E V, ε_{d r}) (\frac{R O R}{R O R + k_{d r}}) (8)$

$G_{1} (B C_{f r}) = \frac{{(B C_{f r})}^{n}}{{(B C_{f r})}^{n} + {(k_{r e v})}^{n}} (9)$

$G_{2} (B C_{f r}) = \frac{(B C_{f r})}{(B C_{f r}) + (k)} (10)$

$F (C_{N}) = \frac{ρ_{c}}{{(1 + k_{c} C_{N})}^{n_{c}}} (11)$

The differential equations are

$\frac{d P_{1}}{d t} = r_{1} F (C_{N}) B C_{f r} - r_{2} P_{1} (12)$

$\frac{d P_{2}}{d t} = r_{2} P_{1} - r_{3} P_{2} (13)$

$\frac{d P_{3}}{d t} = r_{3} P_{2} - r_{4} P_{3} (14)$

$\frac{d P_{4}}{d t} = r_{4} P_{3} + m_{c} P C - d_{4} P_{4} - τ_{c} P_{4} C (15)$

$\frac{d C}{d t} = r_{5} F (C_{N}) B C_{f r} + m_{c} P C - d_{5} C - τ_{c} P_{4} C (16)$

$\frac{d P C}{d t} = τ_{c} P_{4} C - d_{6} P C - m_{c} P C (17)$

$\frac{d P C_{N}}{d t} = r_{6} P C + τ_{N} P_{N} C_{N} - d_{7} P C_{N} - m P C_{N} (18)$

$\frac{d P_{N}}{d t} = - τ_{N} P_{N} C_{N} - d_{P} P_{N} + m P C_{N} (19)$

$\frac{d C_{N}}{d t} = - τ_{N} P_{N} C_{N} - d_{C} C_{N} + m P C_{N} (20)$

$\frac{d B C}{d t} = β_{b c} S - d_{B C} B C (21)$

$\frac{d R E V}{d t} = r_{r e v} G_{1} (B C_{f r}) F (C_{N}) - d_{r e v} R E V (22)$

$\frac{d R O R}{d t} = b_{r o r} + r_{r o r} G_{2} (B C_{f r}) F (C_{N}) - d_{r o r} R O R (23)$

$\frac{d S}{d t} = β + R_{s} (S, R E V) + A_{s} (S, R E V, R O R) - d_{s} S (24)$

$\frac{d T H}{d t} = b_{t h} + R_{t h} (T H, R E V) + A_{t h} (T H, R E V, R O R) - d_{t h} T H (25)$

$\frac{d M A O}{d t} = r_{M} F (C_{N}) B C_{f r} - d_{m} M A O (26)$

$\frac{d D R D 3}{d t} = b_{d r} + R_{d r} (D R D 3, R E V) + A_{d r} (D R D 3, R E V, R O R) - d_{d r} D R D 3 (27)$

The parameters for the Core Circadian Clock are

$\begin{array}{l} K_{d} = 0.02; ρ_{c} = 3; k_{c} = 0.5; n_{c} = 3; r_{1} = 5; r_{2} = 0.45; r_{3} = 0.45; r_{4} = 0.45; \\ m_{c} = 0.5; d_{4} = 0.4; r_{c} 0.5; r_{5} = 5; d_{5} = 0.1; d_{6} = 0.12; r_{6} = 0.75; d_{7} = 0.2; \\ τ_{n} = 0.1; d_{p} = 0.25; d_{c} = 0.2 β_{b c} = 0.1; d_{b c} = 0.1; ξ = 2 / 3 \end{array}$

The parameters for the secondary loop are

$\begin{array}{l} n = 2; k_{r e v} = 0.2; k = 1.5; ρ_{s} = 1; k_{s} = 0.5; ε_{s} = 0; n_{s} = 5.3; α_{s} = 3.7; \\ r_{r e v} = 1.5; d_{r e v} = 0.5; b_{r o r} = 0.1; r_{r o r} = 1.8; d_{r o r} = 0.25; β = 0.9; d_{s} = 3 \end{array}$

The parameters for the equations of the dopamine (DA) elements are

$\begin{array}{l} ρ_{t h} = 1; k_{t h} = 10; ε_{t h} = 0.4; α_{t h} = 1.23; k_{t h} = 1; b_{t h} = 0.85; d_{t h} = 5.6; d_{m} = 0.016 \\ ρ_{d r} = 1; k_{d r} = 10; ε_{d r} = 0.4; n_{d r} = 1; α_{d r} = 0.53; b_{d r} = 0.3; d_{d r} = 3 \end{array}$

Numerical methods used

Bifurcation analysis: Multiple steady-states and oscillatory behavior occur in various situations. Branch and Limit bifurcation points cause multiple steady-states while Hopf bifurcation points produce limit cycles. The MATLAB program MATCONT [41,42] is a commonly used software to locate limit points, branch points, and Hopf bifurcation points. Consider an ODE system

$\dot{x} = f (x, β) (28)$

$x \in R^{n}$ . Defining the matrix A as

$A = [\begin{array}{l} \frac{\partial f_{1}}{\partial x_{1}} \frac{\partial f_{1}}{\partial x_{2}} \frac{\partial f_{1}}{\partial x_{3}} \frac{\partial f_{1}}{\partial x_{4}} .......... \frac{\partial f_{1}}{\partial x_{n}} \frac{\partial f_{1}}{\partial β} \\ \frac{\partial f_{2}}{\partial x_{1}} \frac{\partial f_{2}}{\partial x_{2}} \frac{\partial f_{2}}{\partial x_{3}} \frac{\partial f_{2}}{\partial x_{4}} .......... \frac{\partial f_{2}}{\partial x_{n}} \frac{\partial f_{2}}{\partial β} \\ ........................................................... \\ ........................................................... \\ \frac{\partial f_{n}}{\partial x_{1}} \frac{\partial f_{n}}{\partial x_{2}} \frac{\partial f_{n}}{\partial x_{3}} \frac{\partial f_{n}}{\partial x_{4}} .......... \frac{\partial f_{n}}{\partial x_{n}} \frac{\partial f_{n}}{\partial β} \end{array}] (29)$

β is the bifurcation parameter. The matrix A can be written in a compact form as

$A = [B | \partial f / \partial β] (30)$

The tangent at any point x; ( $v = [v_{1}, v_{2}, v_{3}, v_{4}, .... v_{n + 1}]$ ) must satisfy

Av = 0 (31)

The matrix B must be singular at both limit and branch points. The n+1th component of the tangent vector V_n+1 = 0 at a limit point (LP) and for a branch point (BP) the matrix $[\begin{array}{l} A \\ v^{T} \end{array}]$ must be singular. At a Hopf bifurcation,

$\det (2 f_{x} (x, β) @ I_{n}) = 0 (32)$

@ indicates the bialternate product while I_n is the n-square identity matrix. Hopf bifurcations result in unwanted limit cycles (which in turn cause problems for optimization and control) and should be eliminated. Further details can be found in the works of Kuznetsov [43] and Govaerts [44].

Multiobjective nonlinear model predictive control algorithm

Flores Tlacuahuaz [45] first proposed the Multiobjective nonlinear model predictive control method that does not involve weighting functions, nor does it impose additional constraints on the problem unlike the weighted function or the epsilon correction method [46]. For a set of ODE

$\begin{array}{l} \frac{d x}{d t} = F (x, u) \\ h (x, u) \leq 0 x^{L} \leq x \leq x^{U}; u^{L} \leq u \leq u^{U} (33) \end{array}$

For a final time of tf let p_j (tf) j=1,2,...n be the variables that need to be optimized (maximized or minimized). Simultaneously. n the total number of variables that need to be optimized simultaneously. In this MNLMPC method, dynamic optimization problems that independently minimize/maximize each variable p_j (tf) j=1,2,...n are solved individually. The individual minimization/maximization of each p_j (tf) j=1,2,...n will lead to the values . Then the multiobjective optimal control problem that will be solved is

$\begin{array}{l} \min (\sum_{j = 1}^{n} (p_{j} (t_{f}) - p_{j}^{*}))^{2} \\ s u b j e c t t o \frac{d x}{d t} = F (x, u); h (x, u) \leq 0 (34) \\ x^{L} \leq x \leq x^{U}; u^{L} \leq u \leq u^{U} \end{array}$

This will provide the control values for various times. The first obtained control value is implemented and the rest are ignored. The procedure is repeated until the implemented and the first obtained control values are the same or if the Utopia point p_j (tf) = $p_{j}^{*}$ ; for all j from 1 to n is achieved. The optimization package in Python, Pyomo [47], where the differential equations are automatically converted to algebraic equations will be used. The resulting optimization problem was solved using IPOPT [48]. The obtained solution is confirmed as a global solution with BARON [49]. To summarize the steps of the algorithm are as follows.

Minimize/maximize p_j (t_f) j=1,2,...n. This will lead to the value $p_{j}^{*}$ .
Minimize $(\sum_{j = 1}^{n} (p_{j} (t_{f}) - p_{j}^{*}))^{2}$ . This will provide the control values for various times.
Implement the first obtained control values and discard the remaining.
The steps are repeated until there is an insignificant difference between the implemented and the first obtained value of the control variables or if the Utopia point is achieved.

Results and discussion

The Variables are ordered as

X = [P₁,P₂,P₃,P₄,C,PC,PC_N,PN,C_N,BC,REV,ROR,S,TH,MAO,DRD3]. Bifurcation analysis was performed using the MATLAB program MATCONT. d_s is the bifurcation parameter. The bifurcation analysis revealed the existence of two Hopf bifurcation points given by HA and HB in Figure 1. The variable values at both these bifurcation points are X = (0.668579 0.668579 0.668579 0.369388 2.216329 0.660230 1.197233 1.406020 1.757524 1.230547 0.416188 0.665301 1.230547 0.348979 112.822706 0.575163 1.699851) and X = (0.944861 0.944861 0.944861 0.482939 2.897634 1.128532 2.379765 2.243310 2.804137 3.868766 0.515709 0.723462 3.868766 0.646619 159.445209 1.411344 0.713746), respectively.

The limit cycles that occur at both these points are shown in Figures 2,3. When the bifurcation parameter d_s was modified to $(d_{s} \tanh (d_{s})) / 40;$ both the Hopf bifurcation parameters disappeared. (Figure 4) Sridhar in 2024 [50] explained with several examples how the activation factor involving the tanh activation function where a bifurcation parameter u is replaced by $(u \tanh u / ε)$ ) successfully eliminates the limit cycle causing Hopf bifurcation points. Eliminating the Hopf bifurcation points in the dopamine circadian rhythms model using the tanh factor confirms the results [50].

DRD3 plays an important role in cognition. BC is the protein dimer that drives the transcription of core clock genes. TH (Tyrosine hydroxylase) is the rate-limiting enzyme in dopamine synthesis. MAO (Monoamine oxidase) catalyzes dopamine degradation. Hence, the MNLMPC aims to maximize the final value of the dopamine receptor (DRD3) and BC while minimizing the final values of MAO and TH.

The MNLMPC strategy described previously was used. ds is the control variable. $D R D 3 (t_{f}) + B C (t_{f})$ was first maximized the resulted in a value of 107.1273. Then $T H (t_{f}) + M A O (t_{f})$ was minimized. This resulted in a value of 0.0969. The multiobjective optimal control problem involved the minimization of t ${(D R D 3 (t_{f}) + B C (t_{f}) - 107.1273)}^{2} + {(T H (t_{f}) + M A O (t_{f}) - 0.0969)}^{2}$ . The first obtained control value was implemented and the rest discarded until there was no difference between the implemented and the first obtained control value. This obtained value is referred to as the MNLMPC value.

The MNLMPC value obtained was d_s = 0.0326163. The resulting profiles are shown in Figures 5-7.

Figure 7 shows the d_s versus t profile that demonstrates the existence of spikes. The tanh activation factor effectively eliminates spikes in optimal control profiles [51,52,53]. When the control variable d_s was modified to $(d_{s} \tanh (d_{s})) / 40;$ and the MNLMPC calculations were done again.

The maximization of $D R D 3 (t_{f}) + B C (t_{f})$ produced a value of 107.17. $T H (t_{f}) + M A O (t_{f})$ was then minimized and the value obtained was 0.3567. The multiobjective optimal control problem involved the minimization of ${(D R D 3 (t_{f}) + B C (t_{f}) - 107.17)}^{2} + {(T H (t_{f}) + M A O (t_{f}) - 0.3567)}^{2}$ . The first obtained control value was implemented and the rest discarded until there was no difference between the implemented and the first obtained control value. The MNLMPC value obtained was d_s == 0.6731.

The resulting profiles are shown in Figures 8-10. Figure 10 shows the d_s versus t profile demonstrating that the spikes have been eliminated because of the tanh activation factor.

The multiobjective nonlinear model predictive control is more rigorous than other methods that involve additional unnecessary constraints or equations [46]. The elimination of the Hopf bifurcations avoids the unnecessary limit cycles that may disrupt the circadian rhythm.

Conclusions and future work

The dopamine circadian rhythms model is shown to have two Hopf bifurcations, which cause limit cycles that can disrupt the circadian rhythms. An activation factor involving the tanh function eliminates the limit cycle causing Hopf bifurcations. The high nonlinearity in the dopamine circadian rhythms model also produces spikes when multiobjective nonlinear model predictive calculations are performed. The same activation factor also eliminates the spikes in the control profiles. Future work will involve the performance of a combination of bifurcation analysis and Multiobjective nonlinear model predictive control calculations on more advanced Circadian rhythm models.