Rapid Circadian Entrainment in Models of Circadian Genes Regulation

TL;DR

This paper develops methods for rapid circadian entrainment in biological models by simplifying complex gene regulation equations and applying optimal control techniques, demonstrating effective light input strategies.

Contribution

It introduces a combined approach using PRC, POD, and gradient descent for optimal light control in high-order circadian models, advancing circadian entrainment methods.

Findings

Effective light input strategies for circadian entrainment

Simplification of complex nonlinear models

Impact on gene transcription and protein synthesis

Abstract

The light-based minimum-time circadian entrainment problem for mammals, Neurospora, and Drosophila is studied based on the mathematical models of their circadian gene regulation. These models contain high order nonlinear differential equations. Two model simplification methods are applied to these high-order models: the phase response curves (PRC) and the Principal Orthogonal Decomposition (POD). The variational calculus and a gradient descent algorithm are applied for solving the optimal light input in the high-order models. As the results of the gradient descent algorithm rely heavily on the initial guesses, we use the optimal control of the PRC and the simplified model to initialize the gradient descent algorithm. In this paper, we present: (1) the application of PRC and direct shooting algorithm on high-order nonlinear models; (2) a general process for solving the minimum-time optimal control problem on high-order models; (3) the impacts of minimum-time optimal light on circadian gene transcription and protein synthesis.

Equations

$$\dot{y}=F_{\mathrm{mammal}}(y,lis)=f_0(y)+f_1(y)lis,$$

$$\frac{d[M,F_{\mathrm{N}},F_{\mathrm{C}}{]}^T}{dt}=F_{\mathrm{Neurospora}}(M,F_{\mathrm{N}},F_{\mathrm{C}},v_{\mathrm{s}}),$$

$$\frac{d}{dt}[M_{\mathrm{P}},P_0,P_1,P_2,M_{\mathrm{T}},T_0,T_1,T_2,C,C_{\mathrm{N}}{]}^T$$

$$=F_{\mathrm{Drosophila}}(M_{\mathrm{P}},P_0,P_1,P_2,M_{\mathrm{T}},T_0,T_1,T_2,C,C_{\mathrm{N}},v_{\mathrm{dT}}),$$

$$P_t=P_0+P_1+P_2+C+C_N,T_t=T_0+T_1+T_2+C+C_N.$$

$$\dot{x}=F(x(t),u(t)),$$

$${\phi }_f(t_f,x(t_f))=||x(t_f)-x_{\text{ref}}(t_f)|{|}^2-tol=0,$$

$$u_{\rm ref}(t)=\left\{\begin{array}[]{cl}u_{\rm max},&{\rm mod}(t,24)\in[0,12),\\ u_{\rm min},&{\rm mod}(t,24)\in[12,24).\end{array}\right.$$

$$\dot{x}=F(x(t),u_{\mathrm{ref}}(t)).$$

$$\theta=\left\{\begin{array}[]{cl}{\rm tan}^{-1}(\bar{y}_{2}/\bar{y}_{1}),&{\rm in\ mammalian\ model},\\ {\rm tan}^{-1}(\bar{F}_{\rm C}/\bar{M}),&{\rm in\ Neurospora\ model},\\ {\rm tan}^{-1}(\bar{C}_{\rm N}/\bar{M}_{\rm T}),&{\rm in\ Drosophila\ model},\end{array}\right.$$

$$u_{\rm delay}(\theta,\theta_{\rm ref})=\left\{\begin{array}[]{cl}u_{\rm max},&f(\theta)<0\ {\rm and}\ \Delta\theta\geq\Upsilon,\\ u_{\rm min},&f(\theta)\geq 0\ {\rm and}\ \Delta\theta\geq\Upsilon,\\ u_{\rm ref}(\theta_{\rm ref}),&\Delta\theta<\Upsilon,\end{array}\right.$$

$$u_{\rm advance}(\theta,\theta_{\rm ref})=\left\{\begin{array}[]{cl}u_{\rm max},&f(\theta)>0\ {\rm and}\ \Delta\theta\geq\Upsilon,\\ u_{\rm min},&f(\theta)\leq 0\ {\rm and}\ \Delta\theta\geq\Upsilon,\\ u_{\rm ref}(\theta_{\rm ref}),&\Delta\theta<\Upsilon,\end{array}\right.$$

$$x(t,u(t);x_0)\approx U_{r}z(t,u(t);z_0).$$

$$\dot{z}(t)\approx U_r^{T}F(U_{r}z(t,u(t);z_0),u(t)),\text{where}{z}_0\approx U_r^T{x}_0.$$

$$X_{\mathrm{free}}=\{x_{\mathrm{free}}(t_1)x_{\mathrm{free}}(t_2)...x_{\mathrm{free}}(t_N)\}\in {\mathbb{R}}^{n\times N}.$$

$$X_{\mathrm{free}}=U\Sigma V^{\ast }=[u_1{u}_2...u_n]\Sigma V^{\ast }\approx U_r{\Sigma }_r{V}_r^{\ast }=[u_1{u}_2]{\Sigma }_r{V}_r^{\ast },$$

$$\dot{p}=-\frac{\partial \mathcal{H}}{\partial z}=-\frac{\partial \dot{z}}{\partial z}p.$$

$$u^{\ast }(t)=\arg \underset{u(t)}{\min }\mathcal{H}=\frac{u_{\max }-u_{\min }}{2}[1-\mathrm{sign}(p^T\dot{z})]+u_{\min }.$$

$$p(t_f)=-\frac{\frac{\partial {\phi }_f(t_f,z(t_f))}{\partial z(t_f)}}{\frac{\partial {\phi }_f(t_f,z(t_f))}{\partial t_f}+{\left[\frac{\partial {\phi }_f(t_f,z(t_f))}{\partial z(t_f)}\right]}^T\dot{z}(t_f)}.$$

$$u_{\rm 2nd}(t)=\left\{\begin{array}[]{cl}u_{\rm DSA}(t),&t\leq t_{\rm DSA},\\ u_{\rm ref}(t),&{\rm otherwise}.\end{array}\right.$$

$$J_a(u)={\int }_0^{t_f}1+{\lambda }^T(\tau )[F(x(\tau ),u(\tau ))-\dot{x}(\tau )]d\tau$$

$$+ϵ{\phi }_f(t_f,x(t_f)).$$

$$\dot{\lambda }(t)=-{\lambda }^T\frac{\partial F(x(t),u(t))}{\partial x},\lambda (t_f)=ϵ\frac{\partial {\phi }_f(t_f,x(t_f))}{\partial x},$$

$$ϵ=-\frac{1}{\frac{\partial {\phi }_f(t_f,x(t_f))}{\partial t}+\frac{\partial {\phi }_f(t_f,x(t_f))}{\partial x}F(x(t_f),u(t_f))}.$$

$$\delta J_a{|}_u=\underset{\alpha \to 0}{\lim }\frac{\partial [J_a(u+\alpha \xi )-J_a(u)]}{\partial \alpha }$$

$$={\int }_0^{t_f}{\lambda }^T(\tau )\frac{\partial F(x(\tau ),u(\tau ))}{\partial u}\xi (\tau )d\tau .$$

$${\nabla }_{u(t)}J={\lambda }^T(t)\frac{\partial F(x(t),u(t))}{\partial u}.$$

$$u^{k+1}(t)=\min \left\{\max \left[u^k(t)-{\eta }_u{\nabla }_{u^k(t)}J,u_{\min }\right],u_{\max }\right\},$$

$$\left\{\begin{array}[c]{c}u^{k}(t)=u_{\rm min}\ {\rm when}\ \nabla_{u^{k}(t)}J>0,\\ u^{k}(t)=u_{\rm max}\ {\rm when}\ \nabla_{u^{k}(t)}J<0,\\ u^{k}(t)\in[u_{\rm min},u_{\rm max}]\ {\rm when}\ \nabla_{u^{k}(t)}J=0.\end{array}\right.$$

Full text

I Introduction

††The authors are with the Light Enabled Systems and Application (LESA) Engineering Research Center, Rensselaer Polytechnic Institute, Troy, NY, USA. Research reported in this paper is supported by the NSF through grant number EEC-0812056 and the ARO through grant number W911NF-13-1-0265. The authors would like to thank Dr. Wei Qiao for his contribution to the initial portion of this research as part of his postdoctoral research. Email: {yinj4,juliua2,wenj}@rpi.edu.

A wide range of biological processes in terrestrial species display endogenous oscillations with a period of about 24 hours. We call these processes the circadian rhythm. The circadian rhythm is in close connection with the numerous biological processes, including the sleep-wake cycle, blood pressure, hormone secretion, circadian gene transcription and protein synthesis[1][2].

Maintenance of the circadian rhythm is essential for mammals, insects, and fungi. Being the most powerful influence on the circadian system, light has been used in circadian entrainment in the previous works[3][4]. In view of control theory, the circadian entrainment problem is usually treated as a tracking problem, in which light is used to drive the subject’s circadian state to a reference trajectory[5][6]. The light input for circadian entrainment is usually proposed based on mathematical models of the circadian system. One of the most widely-used circadian models is proposed by Kronauer et al.[7][8]. In this model, the dynamics of the core body temperature is modeled as an oscillator to represent the circadian clock. Light strategies for rapid circadian entrainment on the Kronauer model were given in previous works [9], [10], [11] [12].

The molecular mechanism of the circadian system has also been studied and the circadian genes expression activities are proposed to represent the circadian oscillator. Leloup et al. [13][14] presented a 10th-order ODE Drosophila circadian model, which contains the dynamics of the Period (PER) and Timeless (TIM) proteins. A molecular model of Neurospora was given by Leloup and Gonze [14] in terms of the $frq$ mRNA transcription and FRQ protein synthesis. Numerous models of mammalian circadian rhythm were formulated based on the negative feedback loop in the $Period\ (Per)$ and $Cryptochrome\ (Cry)$ genes transcription[15][16].

In this paper, we seek to solve the minimum-time entrainment problem, i.e., how to use light to synchronize the phase of the circadian rhythm with a reference phase as quickly as possible. We have recently solved this problem for the empirical Kronauer circadian model, which has 3 states (2 oscillatory states and 1 receptor state) [10]. Our approach, based on a functional gradient descent algorithm, searches for the optimal light input for the full-order (3rd order) model using the optimal light inputs of some reduced-order models as initial guesses. We defined two reduced order models: (1) A first-order model that captures only the phase dynamics of the oscillator, and (2) A second-order model that ignores the receptor state dynamics. The first-order model is accurate around the limit cycle of the system, while the second-order model can capture the amplitude dynamics of the oscillation in addition to the phase dynamics. We demonstrated that these initial guesses led to the optimal solution for the full-order system. The circadian models of mammals and Drosophila that are based on circadian genes regulation dynamics have higher orders (7 or more states). For these systems, the optimal light strategy for minimum-time entrainment is hard to calculate directly. Extending the ideas from [10], we propose a general solution procedure listed below:

1. Generate a first-order reduced model that captures the phase dynamics of the oscillation. The first-order dynamics is described by the phase response curve (PRC). Perform PRC-based delay and advance entrainment on the full model, denote the corresponding light inputs as $u_{\rm delay}$ and $u_{\rm advance}$;

2. Use the Principle Orthogonal Decomposition to obtain a second-order oscillator as a reduced model. Perform a direct shooting method [10] to search the optimal light input for the minimum-time entrainment problem, denote it as $u_{\rm DSA}$;

3. Perform a gradient descent algorithm to obtain the locally optimal light strategy $u^{*}$ on the full model. The gradient descent process is initialized by the results in steps 1, 2 and the 12-12-hour daily light-dark cycle $u_{\rm ref}$ in (6).

The circadian models we study in this paper are: A. the 7th-order ODE mammalian model in [17], B. the 3rd-order ODE Neurospora model in[14], C. the 10th-order ODE Drosophila model in[13]. These models are introduced in Section II in details. In Section III, we introduce the theory and implementation of the solution procedure mentioned above. Section IV demonstrates the numerical implementation of entrainment strategies and simulation results. Section V draws some conclusions from the results.

II Mathematic Model and Problem Formulation

II-A Mammalian Model

Based on experimental observations on mice, computational mammalian circadian models were proposed in [15][17] as regulatory loops of $Per$, $Cry$, Bmal1 genes transcription. The model in [17] treats the $Per$ and $Cry$ genes as merely one variable. The 7 variables in this model are: the $Per/Cry$ mRNA ($y_{1}$), cytoplasmic PER/CRY complex ($y_{2}$), nucleus PER/CRY complex ($y_{3}$), cytoplasmic $Bmal1$ mRNA ($y_{4}$), cytoplasmic $Bmal1$ protein ($y_{5}$), nuclear $Bmal1$ ($y_{6}$) and modified nuclear $Bmal1$ ($y_{7}$). The dynamics of the states $y=[y_{1},y_{2},y_{3},y_{4},y_{5},y_{6},y_{7}]^{T}\in\mathbb{R}^{7\times 1}$ is denoted as

[TABLE]

where $lis$ is the light-related input of this model and valued between 0 and a constant value. In terms of [17], we set the maximum value of $lis$ is 0.02 in the presence of light.

II-B Neurospora Model

The circadian rhythm in Neurospora regulates the periodic expression of circadian genes[18]. The circadian model of Neurospora crassa in[14] is proposed based on the negative feedback regulation of the $frq$ gene transcription and FRQ protein synthesis. This model also takes the light effects into account, with the dynamics of this model denoted as

[TABLE]

where $M$ is the $frq$ mRNA, $F_{\rm N}$ and $F_{\rm C}$ are FRQ protein in the nucleus and cytoplasm, respectively. $v_{\rm s}$ denotes the $frq$ transcriptional rate, which is a light-dependent variable, we set $v_{\rm s}=1.6$ in dark and $v_{\rm s}=2$ in full light, i.e., $v_{\rm s}\in[1.6,2]$.

II-C Drosophila Model

Similar to the Neurospora model, the Drosophila circadian model is formulated based on the negative feedback regulation of $per$ and $tim$ genes transcription and PER and TIM proteins synthesis[13][14]. The dynamics of the 10 states of the Drosophila model in [13] are denoted as

[TABLE]

where $M_{\rm P}$ and $M_{\rm T}$ represent the $per$ and $tim$ mRNA, $P_{0},P_{1},P_{2}$ represent three different states of the PER protein, $T_{0},T_{1},T_{2}$ represent three different states of the TIM protein, $C$ and $C_{N}$ represent the PER-TIM complex in the cytoplasm and nucleus, respectively. The total quantities of the PER and TIM proteins are given by:

[TABLE]

The maximum rate of TIM protein degradation $v_{\rm dT}$ is a light-dependent variable. Based on [13], $v_{\rm dT}$ is valued at 2 and 4 in dark and light, respectively.

II-D Problem Formulation

Assume the light-dependent variable in each circadian model is expressed as $u\in\mathbb{R}^{1}$, which is used to represent the light input. In the minimum-time entrainment problem, we want to find the optimal light $u^{*}(t)\in[u_{\rm min},u_{\rm max}]$ to drive the circadian state $x(t)\in\mathbb{R}^{n}$ to a reference trajectory $x_{\rm ref}(t)$ in minimum time given the initial condition $x(0)$ and the system dynamics

[TABLE]

where $n$ is the number of model states. The system dynamics of each model are given as (1), (2) and (3), $u_{\rm min}$ and $u_{\rm max}$ are the upper and lower bounds of the light input. At the final time $t_{f}$, the circadian state should satisfy

[TABLE]

where $tol$ is valued at a small positive scalar and represents the final error between the entraining and reference state.

III Solution Process and Entrainment Strategy

III-A Open-loop Entrainment

The most convenient and widely-used entrainment strategy is applying the daily light-dark cycle on the subject directly. In each model, we assume the reference light $u_{\rm ref}$ is the 12/12 hours light/dark cycle, with a light input given as

[TABLE]

The dynamics of the circadian state in this entrainment process is given as

[TABLE]

The reference trajectory $x_{\rm ref}(t)$ is set as the periodic solution of the differential equation in (7). The light $u(t)=u_{\rm ref}(t)$ is a function of time and does not depend on the circadian state, we call the entrainment by $u_{\rm ref}(t)$ as open-loop entrainment.

III-B PRC and Greedy Delay/Advance Entrainment

PRC is widely used in the study of circadian rhythm [19] and design of light input for circadian entrainment[20]. In this paper, the circadian phase of each model is defined by

[TABLE]

where $\bar{X}$ represents the normalized state $X$, given as $\bar{X}=[2X-(X_{\rm max}+X_{\rm min})]/(X_{\rm max}-X_{\rm min})$, $X_{\rm min}$ and $X_{\rm max}$ are the minimum and maximum values of $X$. The dynamics of the circadian phase is given as $\dot{\theta}=\omega_{0}(\theta)+f(\theta)u$, where $\omega_{0}(\theta)$ is the free running frequency (the frequency in dark). $f(\theta)$ is the phase response curve, which is determined by adding a 30-min light pulse on the free running state (circadian state in dark) and calculating the phase differences with and without light. The results in [20] imply that the optimal control for circadian phase entrainment is either the greedy delay light $u(t)=(u_{\rm max}-u_{\rm min})[1-{\rm sign}(f(\theta))]/2+u_{\rm min}$ or greedy advance light $u(t)=(u_{\rm max}-u_{\rm min})[1+{\rm sign}(f(\theta))]/2+u_{\rm min}$. Note that the amplitude of the circadian state can be quenched or enlarged by the light input. To improve the convergence property, we add a part of reference light as

[TABLE]

[TABLE]

where the phase difference $\Delta\theta$ is defined as $\Delta\theta=\min\{|\theta-\theta_{\rm ref}\pm 2\pi|,|\theta-\theta_{\rm ref}|\}$, and $\Upsilon=1$ rad. Since the reference trajectory keeps synchronized with $u_{\rm ref}(t)$, we represent $u_{\rm ref}(t)$ as $u_{\rm ref}(\theta_{\rm ref})$ in (9) and (10). Assuming the circadian phase of the full model is available during entrainment, we can implement $u_{\rm delay}\ {\rm and}\ u_{\rm advance}$ as feedback control laws.

III-C Model-Order Reduction and Direct Shooting Algorithm

The direct shooting algorithm (DSA) has been introduced in [10] and [11] to solve the minimum-time entrainment problem and search for a globally optimal solution in the second-order Kronauer model. However, this method is impractical for models of third or higher orders. Therefore, before running the shooting algorithm, we reduce the full model into a second-order one. Assume the dynamics of the high-order model is given as (4) with the initial condition $x(0)=x_{0}$. The approximation of state using a right reduced-order base $U_{r}\in\mathbb{R}^{n\times 2}$ is given as

[TABLE]

If $U_{r}^{T}U_{r}=I$, the reduced nonlinear model is given as

[TABLE]

The projection matrix $U_{r}$ of a nonlinear model is usually generated by Principal Orthogonal Decomposition[21], which finds a subspace in $\mathbb{R}^{2}$ that captures information of the original model by minimizing the $\mathcal{L}_{2}$ norm of error $||x(t)-U_{r}z(t)||_{2}$. We simulate the circadian model forward in darkness and generate a free-running trajectory $x_{\rm free}(t)$ in several periods, set a sequence of sampling time $\{t_{1},t_{2},...,t_{N}\}$ and construct a snapshot matrix of the trajectory, denoted as

[TABLE]

By singular value decomposition (SVD), this matrix is decomposed in the following form:

[TABLE]

where $U_{r}\in\mathbb{R}^{n\times 2},\ \Sigma_{r}\in\mathbb{R}^{2\times 2},\ V_{r}\in\mathbb{R}^{N\times 2}$ are truncated SVD matrices which contain the first two principle columns and singular values.

We use the optimal control theory[11] to solve the minimum-time entrainment problem in II-D with a reduced model in (11). The cost function is usually expressed in an integral form as $J=t_{f}=\int_{0}^{t_{f}}1d\tau$. The Hamiltonian of the reduced system is given as $\mathscr{H}=1+p^{T}\dot{z}$, where $p(t)\in\mathbb{R}^{2}$ is the Lagrange multiplier and is usually called the co-state. The dynamics of the co-state is given as

[TABLE]

Based on the Pontryagin Minimum Principle, the optimal control $u^{*}$ satisfies

[TABLE]

The terminal co-state is decided by transversality condition

[TABLE]

The minimum-time entrainment problem is treated as a two-point boundary values problem for solving the optimal initial co-state value $p(0)\in\mathbb{R}^{2}$ given as $p^{*}(0)=\arg\min_{p(0)}t_{f}.$

The solution is subject to the initial condition $z(0)=z_{0}$, the terminal condition $\varphi_{f}(t_{f},z(t_{f}))=0$ and optimality constraints (12), (13), (14). Note that the magnitude of $p(0)=[p_{1}(0),p_{2}(0)]^{T}$ has no effect on the entrainment process and the final result. Thus, we limit our search on a unit circle as $p(0)=[\cos(\phi),\sin(\phi)]^{T}$. The true magnitude of $p(t_{f})$ can later be determined by (14). The direct shooting algorithm only need to search the optimal value of $\phi\in[0,2\pi]$ with a process given in [11]. The light input $u(t)$ corresponding to $p^{*}(0)$ is the optimal light input for minimum-time entrainment of the reduced model, denoted as $u_{\rm DSA}(t)$. We apply $u_{\rm DSA}(t)$ on the full model in the form of (15), where $t_{\rm DSA}$ is the optimal entrainment time of the reduced model.

[TABLE]

III-D Gradient Descent Algorithm

A gradient descent algorithm is used for solving the minimum-time entrainment problem of the full model. The system dynamics and the terminal constraint are plugged into the cost function by introducing multipliers $\lambda(t)\in\mathbb{R}^{n\times 1}$ and $\epsilon\in\mathbb{R}^{1}$[22], the augmented cost function is written as

[TABLE]

For the feasible state $x(t)$ and corresponding time cost $t_{f}$ given $x_{0}$ and $u(t)$, we introduce a small perturbation into $u$, the perturbed light input is given as $u(t)+\alpha\xi(t)$, where $\alpha$ is a scalar. The cost with variation in the light input is given as $J_{a}(u+\alpha\xi)$. Note that $\lambda(t)$ and $\epsilon$ have no effects on the cost value. For simplification, $\lambda(t)$ and $\epsilon$ are defined as:

[TABLE]

The first order variation of cost with respect to light is

[TABLE]

As $J_{a}=J$, the gradient descent direction of the light is

[TABLE]

Steps of the gradient descent process for solving the optimal light input of the minimum-time circadian entrainment problem are listed below:

Lemma: The gradient descent process reaches a fixed point $u^{*}(t)=u^{k+1}(t)=u^{k}(t)$ for any positive step $\eta_{u}$ if and only if the following conditions are satisfied for $\forall t\in[0,t_{f}^{k}]$:

[TABLE]

The result in Algorithm 8 relies heavily on the initial guess. To get a locally optimal solution with small entrainment time, we choose four initial guesses for GDA: a) $u^{0}(t)=u_{\rm ref}(t)$, b) $u^{0}(t)=u_{\rm advance}(t)$, c) $u^{0}(t)=u_{\rm delay}(t)$, d) $u^{0}(t)=u_{\rm 2nd}(t)$.

IV Numerical Implementation

Consider the case that the subject is entrained under the reference light at the beginning. A sudden time shift $\Delta_{\rm shift}$ in the external light-dark cycle occurs at $t=0$ and the subject loses synchronization with the local environment, that is, $x(0)=x_{\rm ref}(\Delta_{\rm shift})\neq x_{\rm ref}(0)$, $\Delta_{\rm shift}\in\{1,2,...,23\}$. We want to re-entrain the subject as quickly as possible.

IV-A Mammalian Model

For the mammalian model, the circadian state $x(t)=y(t)$ and the light input $u=lis$. Fig. 2 shows the gradient descent results with four initial guesses (in Fig. 2a) and various entrainment time costs (in Fig. 2b). The PRC-based delay-advance threshold occurs at about 13 hours shift, the delay/advance entrainment works better than $u_{\rm ref}$ and $u_{\rm 2nd}$ in the cases with a time shift of 8-16 hours. The GDA results show that different initial lights may converge to different locally optimal solutions. In the 11 hours shift case, the GDA initialized by $u_{\rm advance}(t)$ reduces the time cost from 260 to 154 hours, while that initialized by $u_{\rm delay}(t)$ decreases the cost from 210 to 164 hours.

Fig. 3 shows 8, 12 and 16-hour shifts cases by the optimal light and time evolution of the $Per/Cry$ mRNA and protein during the optimal entrainment process. The lower sub-figures in Fig. 3(a), 3(c), 3(e) compare the gradient $\nabla_{lis}J$ with the optimal light $lis^{*}(t)$, implying that the light $lis^{*}(t)$ in these figures follows the optimality condition mentioned in (22). Note that the optimal controls in these cases are all bang-bang control, i.e., $lis^{*}(t)$ is either maximum or minimum. Fig. 3(b), 3(d) and 3(f) show the $Per/Cry$ mRNA–PER/CRY protein cycle during the optimal entrainment process. Define the time shift corresponding to the maximum time cost in the optimal entrainment as $\Delta_{\rm threshold}$. In this mammalian model, $\Delta_{\rm threshold}\approx 11$ hours, we can seen that when $\Delta_{\rm shift}>\Delta_{\rm threshold}$, the amplitude of $Per$-$Cry$ mRNA–protein oscillator is quenched during entrainment, i.e., for 12 and 16 hours shifts; while in the case with 8 hours shift, the amplitude of mRNA-protein cycle is enlarged, the optimal light tend to enhance the $Per/Cry$ mRNA transcription and PER/CRY protein synthesis in this case.

IV-B Neurospora Model

In the Neurospora model, $x=[M,F_{\rm C},F_{\rm N}]^{T}$ and $u(t)=v_{\rm s}(t)$. Fig. 4(a) shows the GDA results in the Neurospora model, we can observe that the DSA results from the reduced models are the best initial guess of the GDA as the GDA initialized by $u_{\rm 2nd}$ always reaches the best solution among all four initial guesses in every case. Fig. 4(b) demonstrates the entrainment time of four strategies in the Neurospora model. Compared with the open-loop entrainment and entrainment by $u_{\rm 2nd}$, the delay/advance strategy and optimal light dramatically decrease the time cost, especially in cases with $\Delta_{\rm shift}\in[3,20]$ hours. The maximum time cost in the optimal entrainment occurs at $\Delta_{\rm threshold}=$4 hours shift. This implies the Neurospora model is an asymmetrical one in which advancing 18 hours is faster than delaying 6 hours.

Fig. 5 plots 8, 12 and 16-hour shifts optimal entrainment cases of the Neurospora model. Similar to the results in the mammalian model, the optimal control in the Neurospora model is a bang-bang control. From the three cases ($\Delta_{\rm shift}>\Delta_{\rm threshold}$) shown in Fig. 5, we observe that during the optimal entrainment process, the amplitudes of the Neurospora oscillator are all quenched, the optimal light decreases the rate of $frq$ mRNA transcription and FRQ protein synthesis.

IV-C Drosophila Model

Fig. 6(a) shows that $u_{\rm 2nd}$ is the best initial guess for GDA in the Drosophila model. The entrainment time costs in Fig. 6(b) show the maximum time cost of the optimal entrainment occurs at $\Delta_{\rm threshold}=$14 hours shift, with a value of 50 hours. In Fig. 7, we use the $C_{\rm N}$-$M_{\rm T}$-$T_{t}$ cycle to demonstrate the circadian oscillator of Drosophila. Similar with the mammalian and Neurospora model, the amplitude of circadian oscillator is enlarged when $\Delta_{\rm shift}<\Delta_{\rm threshold}$ and quenched when $\Delta_{\rm shift}>\Delta_{\rm threshold}$ during the minimum-time optimal entrainment. These phenomena imply that, in these models, the minimum-time optimal light delays and advances the circadian rhythm by enhancing and inhibiting the clock gene transcription and protein synthesis, respectively.

Note that in the 8-hour shift case, the GDA converge to a light solution that is no longer a bang-off one with $\nabla_{v_{\rm dT}(t)}J=0$ in some intervals. We call this case as a singular case.

V Discussions and Conclusions

In this paper, we solve the minimum-time entrainment problem in the high order mammalian, Neurospora, Drosophila circadian genes regulation model. We first obtain the optimal lighting strategies for the 1st- and 2nd-order reduced models. The gradient descent algorithm then searches for (locally) optimal solutions with those from the reduced models as initial guesses.

Impacts of initial guesses on gradient descent results: The optimal light from the 2nd-order model works better on the GDA than other initial guesses in the Drosophila and Neurospora model. However, the results in the mammalian model show that the greedy advance strategy is the best initial guess for the GDA in a certain case. Therefore, in the solution procedure, we should initialize the GDA by both the greedy delay/advance strategy and the 2nd-order light strategy, and choose the one with the minimum time cost.

Impacts of the optimal light on circadian gene-protein: Among all three models in this paper, the amplitude of the gene-protein oscillator is quenched or enlarged in every case. The results demonstrate that, during optimal entrainment, circadian delay and advance are closely linked with enlargement and quenching of the gene-protein oscillator in these models.

Compared with the minimum-time entrainment, the delay/advance strategy takes longer time. However, the entrainment time of this strategy is less than $u_{\rm ref}$ and $u_{\rm 2nd}$ in mammalian and Neurospora models. This strategy is attractive as a feedback controller if the measurement of the circadian phase is available during the entrainment process.