Bifurcations and control studies in Circadian Rhythms in <em>Drosophila</em>

Lakshmi N Sridhar*

Review Article

More Information

Submitted: January 17, 2025 | Approved: February 06, 2025 | Published: February 07, 2025

How to cite this article: Sridhar LN. Bifurcations and control studies in Circadian Rhythms in Drosophila. Arch Case Rep. 2025; 9(2): 054-065. Available from:
https://dx.doi.org/10.29328/journal.acr.1001128

DOI: 10.29328/journal.acr.1001128

Copyright license: © 2025 Sridhar LN. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Keywords: Bifurcation; Optimization; Control; Drosophila

Bifurcations and control studies in Circadian Rhythms in Drosophila

Lakshmi N Sridhar*

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

*Address for Correspondence: Lakshmi N Sridhar, Chemical Engineering Department, University of Puerto Rico, Mayaguez, PR 00681, Puerto Rico, Email: [email protected]

Abstract

Bifurcation analysis and Multiobjective Nonlinear Model Predictive Control (MNLMPC) calculations were performed on a model of circadian oscillations of the period (PER) and timeless (TIM) proteins in Drosophila. The MATLAB program MATCONT was used to perform the bifurcation analysis. The optimization language PYOMO was used along with the state-of-the-art global optimization solvers IPOPT and BARON for the MNLMPC calculations. The bifurcation analysis revealed oscillation causing Hopf bifurcations while the MNLMPC calculations revealed the existence of spikes in the control profiles. Both Hopf bifurcation points and the control profile spikes were eliminated using an activation factor involving the hyperbolic tangent function.

Background

Decroly, et al. [1] discovered the chaotic behavior in multiply regulated biochemical systems. Alamgir and Epstein [2] studied the oscillations in coupled chemical oscillators of the chlorite-bromate-iodide system. Aronson, et al. [3] studied the negative feedback defining a circadian clock. Baylies, et al. [4] conducted genetic, molecular, and cellular studies of the per locus and its products in Drosophila melanogaster. Crosthwaite, et al. [5] investigated the photo responses and the origins of circadian rhythmicity. Curtin, et al. [6] demonstrated how the temporally regulated nuclear entry of the Drosophila period protein contributes to the circadian clock. Dunlap [7] performed a genetic and molecular analysis of circadian rhythms. Edery, et al. [8.9] studied the phase shifting of the circadian clock by induction of the Drosophila period protein and its Temporal phosphorylation. Edmunds [10] discussed models and mechanisms for Circadian Timekeeping.

Eskin, et al. [11] investigated the requirement for protein synthesis to regulate the circadian rhythm by melatonin. Gekakis, et al. [12] demonstrated the Defective interaction between timeless protein and long-period mutant. Goldbeter [13] developed a model for circadian oscillations in the Drosophila period (PER) protein. Goldbeter [14] wrote a textbook in Biochemical Oscillations and Cellular Rhythms. Goodwin [15] minvestigated the oscillatory behavior in enzymatic control processes. Hall and co-workers [16-18] studied the molecular effects on biological rhythms. Hong and Tyson [19] investigated temperature compensation of the circadian rhythm in Drosophila based on dimerization of the PER protein. Huang, et al. [20] studied PER protein interactions and temperature compensation of a circadian clock in Drosophila. Khalsa, et al. [21] studied techniques for stopping the circadian pacemaker with inhibitors of protein synthesis Konopka and Benzer [22] investigated clock mutants of Drosophila melanogaster. Lee, et al. [23] discuss strategies for resetting the Drosophila clock by photic regulation of PER. and a PER-TIM complex. Leloup and Goldbeter [24] researched temperature compensation of circadian rhythms. Marus, et al. [25] studied the effect of constant light and circadian entrainment of perS flies. Myers, et al. [26] studied the light-induced degradation of timeless and entrainment of the Drosophila circadian clock. Qiu and Hardin [27] showed that the per mRNA cycling is locked to lights-off under photoperiodic conditions that support the circadian feedback loop function. Rosbash [28] studied the molecular control of circadian rhythms. Rutila, et al. [29] discuss t the tim^sl mutant of the Drosophila rhythm gene timeless manifests allele-specific interactions with period gene mutants. Saez and Young [30] discuss the regulation of nuclear entry of the Drosophila clock proteins period and timeless. Saunders, et al. [31] discuss the Light-pulse phase response curves for the locomotor activity rhythm in period mutants of Drosophila melanogaster. Sehgal, et al. [32] studied the loss of circadian behavioral rhythms and per RNA oscillations in the Drosophila mutant timeless. So and Rosbash [33] showed that post-transcriptional regulation contributes to Drosophila clock gene mRNA cycling. Taylor, et al. [34] demonstrate that inhibitors of protein synthesis on 80S ribosomes phase shift the Gonyaulax clock. Vosshall, et al. [35] show the existence of a block in nuclear localization of period protein by a second clock mutation, timeless. Young, et al. [36] discuss the molecular anatomy of a light-sensitive circadian pacemaker in Drosophila.

Zeng and co-workers [37,38] conducted further research about the Drosophila circadian clock. Reitz, et al. [39], Mia, et al. [40], Parcha, et al. [41] and Teeple, et al. [42] studied the relationships between obesity and Circadian rhythms, while Tobeiha, et al. [43] and Festus, et al. [44] studied the effects of cardiovascular health on Circadian rhythms.

Leloup JC, Goldbeter [45] developed a model for circadian rhythms in Drosophila incorporating the formation of a complex between the PER and TIM proteins. In this work a) bifurcation analysis is performed on the Leloup Goldeter model for circadian rhythms in Drosophila to identify and eliminate the Hopf bifurcations and b) Multiobjective NonlinearMmodel Predictive Control(MNLMPC) calculations are performed to ensure maximum production of the required product while simultaneously minimizing all the unwanted by-products.

Motivation and objectives

Circadian rhythms are endogenous limit cycle oscillations characterized for 24 hours. They constitute the biological rhythms with the longest period known to be generated at the molecular level. Oscillations are caused by Hopf bifurcation point. Limit cycle oscillations occur because of Hopf bifurcations. These bifurcations cause spikes in control profiles when one tries to control the process to achieve maximum benefit. These spikes hinder optimization and control tasks. The motivation of this work is to eliminate the limit cycle causing Hopf bifurcations and the spikes in control profiles when dynamic optimization tasks are performed. Hence, the main objectives of this work are a) Identify the Hopf bifurcation points computationally perform bifurcation analysis using MATCONT (a Matlab software for performing bifurcation analysis) b) Use an activation factor involving the tanh function to eliminate these Hopf bifurcation points c) Perform Multiobjective Nonlinear Model Predictive Control(MNLMPC) calculations for the Drosophila circadian rhythms model (Leloup Goldeter; [45]) d) Demonstrate the existence of spikes in the control profile when these calculations are performed) Use the tanh function’s activation factor to eliminate the control profile spikes. This paper is organized as follows. First, the Drosophila circadian rhythms model (Leloup Goldeter; [45]) details are presented. This is followed by describing the bifurcation analysis and Multiobjective Nonlinear Model Predictive Control(MNLMPC) procedures. The results and discussion section is then presented followed by the conclusions.

Drosophila circadian rhythms model

In this model, PER stands for period protein TIM stands for timeless protein. M_p, P₀, P₁, P₂ stand for cytosolic concentration of PER mRNA, unphosphorylated PER, monophosphorylated PER and bisphosphorylated PER. M_T, T₀, T₁, T₂ stand for cytosolic concentration of TIM mRNA, unphosphorylated TIM, monophosphorylated TIM and bisphosphorylated. C represents the complex formed by the degradation of P₂ and T₂ while CN represents the nuclear form of the PER-TIM complex. The mRNAs are synthesized at rates of v_SP, v_ST with rate constants K_sp, K_sT. The subsequent degradation takes place at rates v_mP, v_mT with rate constants K_mp, K_mT. k_Ip, k_2p, k_3p, k_4p are the kinase rate constants while k_1T, k_2T, k_3T, k_4T are the phosphatase rate constants. v_1p, v_2p, v_3p, v_4p represents the maximum kinase rate while v_1T, v_2T, v_3T, v_4T represents the maximum phosphatase rate. More details can be found in Leloup and Goldeter [45]. The base parameter values are

$\begin{array}{l} v_{s P} = v_{s T} = 1 n M h^{- 1}; v_{m P} = v_{m T} = 0.7 n M h^{- 1}; v_{m P} = v_{m T} = 0.7 n M h^{- 1} \\ k s P = k s T = 0.9 h^{- 1}; v_{d P} = v_{d T} = 2 n M h^{- 1}; k 1 = 0.6 h - 1, k 2 = 0.2 h - 1, k 3 = 1 . 2 n M^{- 1} h^{- 1} \\ k_{4} = 0.6 h^{- 1}, K_{I P} = K_{I T} = 1 n M, K_{d P} = K_{d T} = 0.2 n M, n = 4 \\ K_{1 P} = K_{1 T} = K_{2 P} = K_{2 T} = K_{3 P} = K_{3 T} = K_{4 P} = K_{4 T} = 2 n M, k_{d} = k_{d C} = k d N = 0.01 h^{- 1} \\ V_{1 P} = V_{1 T} = 8 n M h^{- 1}, V_{2 P} = V_{2 T} = 1 n M h^{- 1}, V_{3} P = V_{3 T} = 8 n M h^{- 1}, V 4 P = V 4 T = 1 n M h^{- 1} \end{array}$

The equations that constitute this model are

$\frac{d M_{P}}{d t} = V_{S P} (\frac{K_{I P}^{n}}{K_{I P}^{n} + C_{N}^{n}}) - V_{m P} (\frac{M_{P}}{K_{M P} + M_{P}}) - k_{d} M_{P}$ (1)

$\frac{d P_{0}}{d t} = k_{s P} M_{P} - V_{1 P} (\frac{P_{0}}{K_{1 P} + P_{0}}) + V_{2 P} (\frac{P_{1}}{K_{1 P} + P_{1}}) - k_{d} P_{0}$ (2)

$\frac{d P_{1}}{d t} = V_{1 P} (\frac{P_{0}}{K_{1 P} + P_{0}}) - V_{2 P} (\frac{P_{1}}{K_{2 P} + P_{1}}) - V_{3 P} (\frac{P_{1}}{K_{3 P} + P_{1}}) + V_{4 P} (\frac{P_{2}}{K_{4 P} + P_{2}}) - k_{d} P_{1}$ (3)

$\frac{d T_{2}}{d t} = V_{3 T} (\frac{T_{1}}{K_{3 T} + T_{1}}) - V_{4 T} (\frac{T_{2}}{K_{4 T} + T_{2}}) - V_{d T} (\frac{T_{2}}{K_{d T} + T_{2}}) + k_{4} C - k_{3} P_{2} T_{2} - k_{d} T_{2}$ (4)

$\frac{d P_{2}}{d t} = V_{3 P} (\frac{P_{1}}{K_{3 P} + P_{1}}) - V_{4 P} (\frac{P_{2}}{K_{4 P} + P_{2}}) - V_{d P} (\frac{P_{2}}{K_{d P} + P_{2}}) + k_{4} C - k_{3} P_{2} T_{2} - k_{d} P_{2}$ (5)

$\frac{d M_{T}}{d t} = V_{S T} (\frac{K_{I T}^{n}}{K_{I T}^{n} + C_{N}^{n}}) - V_{m T} (\frac{M_{T}}{K_{M T} + M_{T}}) - k_{d} M_{T}$ (6)

$\frac{d T_{0}}{d t} = k_{s T} M_{T} - V_{1 T} (\frac{T_{0}}{K_{1 T} + T_{0}}) + V_{2 T} (\frac{T_{1}}{K_{1 T} + T_{1}}) - k_{d} T_{0}$ (7)

$\frac{d T_{1}}{d t} = V_{1 T} (\frac{T_{0}}{K_{1 T} + T_{0}}) - V_{2 T} (\frac{T_{1}}{K_{2 T} + T_{1}}) - V_{3 T} (\frac{T_{1}}{K_{3 T} + T_{1}}) + V_{4 T} (\frac{T_{2}}{K_{4 T} + T_{2}}) - k_{d} T_{1}$ (8)

$\frac{d C}{d t} = - k_{4} C + k_{3} P_{2} T_{2} + k_{2} C_{N} - k_{1} C - k_{d C} C$ (9)

$\frac{d C_{N}}{d t} =_{2} - k_{2} C_{N} + k_{1} C - k_{d N} C_{N}$ (10)

Bifurcation analysis

Multiple steady-states and oscillatory behavior occur in various situations. Multiple steady states occur because of t Branch and Limit bifurcation points cause multiple steady-states. Hopf bifurcation points produce oscillatory behavior. Ions and limit cycles. The MATLAB program MATCONT. (Dhooge Govearts, and Kuznetsov [46], Dhooge Govearts, Kuznetsov, Mestrom and Riet, [47]) is commonly used software to locate limit points, branch points, and Hopf bifurcation points. Consider an ODE system

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

x ∈ Rn . 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}]$ (12)

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

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

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

Av = 0 (14)

The matrix B must be singular at both limit and branch points. The n+1th component of the tangent vector Vn+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$ (15)

@ indicates the bialternate product while In is the n-square identity matrix. Hopf bifurcations cause unwanted oscillatory behavior and should be eliminated because oscillations make optimization and control tasks very difficult. More details can be found in Kuznetsov [48,49] and Govaerts [50].

Multiobjective nonlinear model predictive control

Flores Tlacuahuaz [51] 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(Miettinen, [52]). For a 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} \end{array}$ (16)

let $\sum_{t_{i = 0}}^{t_{i} = t_{f}} p_{j} (t_{i})$ (j = 12..n); be the variables that need to be minimized/maximized simultaneously, tf being the final time value, and 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 $\sum_{t_{i = 0}}^{t_{i} = t_{f}} p_{j} (t_{i})$ are solved individually. The minimization/maximization of each $\sum_{t_{i = 0}}^{t_{i} = t_{f}} p_{j} (t_{i})$ will lead to the values $p_{j}^{*}$ . Then the optimization problem that will be solved is

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

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 ( $\sum_{t_{i = 0}}^{t_{i} = t_{f}} p_{j} (t_{i}) = p_{j}^{*}$ ; for all j) is achieved. The optimization package in Python, Pyomo (Hart, et al. [53]), where the differential equations are automatically converted to algebraic equations will be used. The resulting optimization problem was solved using IPOPT (Wächter and Biegler [54]). The obtained solution is confirmed as a global solution with BARON (Tawarmalani, M. and N. V. Sahinidis [55]).

Results and discussion

For the three bifurcation parameters kSp, kdN, vmP; the bifurcation analysis revealed Hopf bifurcation points that disappeared when the bifurcation parameters were modified with the activation factor involving the tanh function. The bifurcation parameters kSp, kdN, each resulted in one Hopf bifurcation point, while vmP resulted in two Hopf bifurcation points. The bifurcation parameters were treated as time dependent variables while the other parameters were the base parameters. Whenks ksP was the bifurcation parameter the Hopf bifurcation point was found at the point x = (Mp, P0, P1, P2, MT, T0, T1, T2, C, CN, ksP) = (0.311096 0.847336 1.070162 6.587611 0.311096 0.082204 0.076131 0.030362 0.375862 1.073892 6.559451). The limit cycle resulting from this Hopf bifurcation point is shown in Figure 1a. The Hopf bifurcation point and the limit cycle disappeared when ksP was modified to ksP tanh(ksP) / 10 (Figure 1b).

ISSN: 2637-3793