<title>Introduction</title> <p>Control theoretic tools have been used to model mRNA transcriptional/translational regulatory feedback mechanisms ##UREF##0##[1]##, to analyze nonlinear phenomena ##REF##17158515##[2]##,##REF##15340155##[3]##, and to control complex biological behavior ##UREF##1##[4]##,##UREF##2##[5]##. In our research, we couple systems theoretic tools (such as sensitivity analysis) with model predictive control, to better address phase resetting properties of nonlinear biological oscillators. Our work aims to alleviate circadian-related disorders (such as jet lag and advanced/delayed sleep phase syndromes) by investigating the phase resetting properties of an example circadian mathematical model. More specifically, we manipulate multiple control inputs (or target parameters) to drive the dynamic behavior of the system.</p> <p>Many researchers have shown that the systematic application of light pulses may reset the phase of circadian clocks. This light pulse (input) to induced phase-shift (output) mapping is most notably characterized by the phase response curve (PRC). Daan and Pittendrigh studied the PRC to establish a relationship among circadian behavior (nocturnal vs. diurnal activity), free-running period, and maximum phase advance/delay ##UREF##3##[6]##. The free-running period of an organism reflects its circadian behavior without the influence of entrainment factors such as environmental light∶dark cycles. The free-running period of nocturnal animals, for instance, is often less than 24 hours such that dusk triggers a phase delay and the onset of activity. Conversely, diurnal animals often exhibit free-running periods greater than 24 hours such that dawn triggers a phase advance and the onset of activity ##UREF##3##[6]##. Other researchers have made use of PRCs to establish light as a means to accelerate circadian entrainment ##REF##12398256##[7]##, or as a means to start, stop, and reset the phase of simplified circadian models ##REF##15363673##[8]##–##REF##17952058##[11]##.</p> <p>In a previous study, we develop a closed-loop nonlinear model predictive control (MPC) algorithm that minimizes the phase difference between a reference and a controlled system (each modeled as a single deterministic oscillator) through the systematic application of continuous light. Through use of MPC, circadian phase is recovered in almost half the time required by the natural open-loop sun cycles ##UREF##1##[4]##. Next, we investigated how the MPC algorithm's tuning parameters might affect the model's phase resetting dynamics ##UREF##5##[12]##. Here, we make use of sensitivity analysis to identify additional control inputs (or drug targets) that, when used by the MPC algorithm, outperform light-based circadian phase resetting. The target identification of single and multiple control inputs, coupled with the analysis of their respective performance, parallels efforts in the pharmaceutical industry to yield the greatest behavioral response with respect to the smallest system perturbation. In other words, our methodology may be used to identify optimal (and arguably non-intuitive) drug targets for therapy.</p> <p>To establish an upper bound relating to the time required to recover phase differences, we begin by evaluating the open-loop control algorithm in the <italic>Open-Loop Phase Recovery</italic> section. The identification and manipulation of a set of single, dual, and triple control inputs are then used to minimize phase recovery dynamics of a wild-type circadian system (as described in the <italic>Single</italic>, <italic>Dual</italic>, and <italic>Triple Target Phase Resetting</italic> sections, respectively). This case is most similar to resetting a healthy organism's phase when subject to an environmental disturbance such as jet lag. In the <italic>Short and Long Period Mutants</italic> section, we further investigate how MPC may be used to alleviate chronic circadian disorders. More specifically, we apply the algorithm to circadian oscillator models that exhibit either short or long-period mutant phenotypes. Results suggest that organisms with such syndromes may track regular 24 hour rhythms through the systematic application of light. Our findings support this unique application of systematic drug target identification coupled with model predictive control for use in medicine and pharmacology (see the <xref ref-type="sec" rid="s3">
<italic>Discussion</italic>
</xref> section). In the <xref ref-type="sec" rid="s4">
<italic>Methods</italic>
</xref> section, we describe the employed model predictive control algorithm and the state-based sensitivity analysis used to identify single and multiple parametric control inputs.</p>
<title>Methods</title> <p>A 10-state <italic>Drosophila melanogaster</italic> circadian limit cycle oscillator serves as the model system. This model consists of two coupled auto-regulatory transcription/translation negative feedback loops that characterize <italic>period</italic> and <italic>timeless</italic> gene and protein dynamics ##REF##9486845##[13]##. As demonstrated in previous work, the MPC algorithm may be applied to any stable limit cycle oscillator, including a more complex <italic>Mus musculus</italic> model ##UREF##1##[4]##. Thus, we describe the example system as a general set of nonlinear ordinary differential equations with time <italic>t</italic>, <italic>n</italic>-length state vector <bold>x</bold>(<italic>t</italic>), environmental light input <italic>l</italic>(<italic>t</italic>), additional control inputs <bold>u</bold>(<italic>t</italic>), and system dynamics <bold>f</bold>(<bold>x</bold>(<italic>t</italic>), <italic>l</italic>(<italic>t</italic>), <bold>u</bold>(<italic>t</italic>)):Given that both environmental light and additional control variables may be modeled as multiplicative inputs, the nominal wild-type (sun-cycle entrained) case requires <italic>u</italic>(<italic>t</italic>) = 1, while <italic>l</italic>(<italic>t</italic>) oscillates as a square wave with a frequency of 24 hours, between values 1 and 2. For consistency, the natural sun-cycle environment (or reference) is characterized by the nominal <italic>Drosophila melanogaster</italic> model and denoted by <bold>r</bold>(<italic>t</italic>). This reference is pre-entrained to normal 24 hour light∶dark cycles and is not subject to additional control inputs.</p> <title>Model Predictive Control</title> <p>Model predictive control ##UREF##11##[38]## is used to increase the re-synchronization or entrainment rate of circadian oscillators through the systematic application of specified control inputs. The algorithm follows a sample and hold strategy, updating the prediction and control input every <italic>t<sub>s</sub></italic> = 2 hours, where the discrete time index , such that a function <italic>g</italic>(<italic>kt<sub>s</sub></italic>) = <italic>g</italic>[<italic>k</italic>]. For simplicity, we refer to <italic>k</italic> as being equivalent to  and ignore its rounding component. The manipulated control variable, <italic>u</italic>[<italic>k</italic>], optimizes an open-loop performance objective on a time interval extending from the current time  to the current time plus a prediction horizon of <italic>P</italic> = 48 hours, where . This horizon allows the algorithm to take control action at the current time in response to a forecasted error. The move horizon, <italic>M</italic> = 12 hours, limits the number of control inputs within the prediction horizon such that <italic>u</italic>[<italic>k</italic>] spans a time interval . Beyond  hours of simulation, the predictive model defaults to <italic>u</italic>[<italic>k</italic>] = 1. Future behaviors for a variety of control inputs are computed according to the mathematical model of the system ##REF##9486845##[13]##.</p> <p>The efficacy of the algorithm was evaluated with respect to a sample and hold time interval of 1, 2, and 3 hours (reflecting a move horizon of 3, 6, and 9 hours, respectively). Although shorter light pulses offer a more dynamic manipulated variable profile, it shortens the move horizon and may reduce the utility of model predictive control. Conversely, a longer pulse may reduce the possible control profiles since extended exposure to light leads to arrhythmic behavior ##REF##15746913##[39]##. Thus, we set the sampling rate to 2 hours.</p> <p>The fitness function penalizes the normalized predicted state error between the reference and controlled trajectories, <bold>ē</bold>[<italic>k</italic>], and the net control, <bold>ū</bold>[<italic>k</italic>], over the prediction horizon. The system output used to evaluate circadian performance (or, phase entrainment) is the trajectory defined by the total <italic>period</italic> and <italic>timeless</italic> protein complex concentrations. This state error, <bold>e</bold>[<italic>k</italic>], is normalized with respect to the nominal amplitude of oscillation while the time dependent control input, <bold>u</bold>[<italic>k</italic>], is normalized with respect to the nominal set of values, 1:where the state dynamics <bold>r</bold>[<italic>k</italic>] characterize the nominal reference. Note that the vector <bold>e</bold>[<italic>k</italic>] is , while the matrix <bold>ū</bold>[<italic>k</italic>] is <italic>m</italic>×<italic>c</italic> ( and <italic>c</italic> denotes the number of control inputs).</p> <p>To avoid penalizing transient effects, the state error is weighted uniformly over the move horizon (reflected in the first <italic>m</italic> diagonal values of the <italic>p</italic>×<italic>p</italic> matrix <bold>Q</bold>), and with increasing weight of slope 2 over the prediction horizon (reflected in the <italic>p</italic>−<italic>m</italic> to <italic>p</italic> diagonal values of <bold>Q</bold>). The cost of applying a light input is weighted uniformly with a magnitude of 100 as reflected in the diagonal values of the <italic>m</italic>×<italic>m</italic> matrix <bold>R</bold>. We can afford to be conservative with the cost of control in the wild-type case, since we can ensure that the lack of control (the open-loop algorithm) will eventually entrain the system. The values contained in <bold>R</bold> will be re-evaluated when the algorithm is designed to entrain mutant phenotype models. The performance of an <italic>m</italic>-length control input is measured byOnly the first move of the lowest cost control sequence evaluated at time <italic>k</italic>, , is implemented. Therefore, the sequence of actually implemented control moves may differ significantly from the sequence of control moves calculated at a particular time step. This discrepancy disappears as the prediction and move horizons near infinity. Feedback is incorporated by using the next measurement to update the optimization problem. Once the controlled state trajectories converge to within 15% of the reference state trajectories, the system is considered to have recovered its phase in <italic>T<sub>r</sub></italic> = min<italic><sub>k</sub></italic>[|<italic>e</italic>[<italic>k</italic>]|<sub>∞</sub>≤0.15] hours. At this point, the algorithm defaults to no control since nominal light∶dark cycles will keep the system synchronized to the new environment. Optimization of the phase synchronizing control sequences is completed through use of a genetic algorithm ##UREF##12##[40]##–##UREF##14##[42]##.</p> <title>Sensitivity Analysis</title> <p>Parametric sensitivity analysis quantifies the relative change of system behavior with respect to an isolated parametric perturbation. Parametric <italic>state</italic> sensitivity analysis assigns a value to each system parameter that defines how its perturbation affects state dynamics: . This tool is often used to identify the robustness and fragility tradeoffs of regulatory structures ##REF##15340155##[3]##, and may be tailored to evaluate specific output performance such as period, amplitude, or phase characteristics ##REF##17158515##[2]##.</p> <p>Assuming the model has <italic>n</italic> states and <italic>ρ</italic> parameters, the FIM is a <italic>ρ</italic>×<italic>ρ</italic> matrix describing how any two parametric perturbations might affect state dynamics. More notably, the diagonal values of the FIM describe how any single parameter may affect state dynamics. As a result, we sort the values of the FIM from greatest to least magnitude and choose the top three individual parameters (reflected by the sorted diagonal values) and top three pairs of parameters (reflected by the sorted off-diagonals) whose perturbations yield the greatest change in output.</p> <p>We further analyze the FIM via the singular value decomposition ##UREF##15##[43]##. Assuming FIM = <italic>F</italic>, it may be decomposed as <italic>F</italic> = <italic>UΣV<sup>T</sup></italic>, where Σ is an <italic>n</italic> by <italic>p</italic> diagonal matrix of non-negative singular values, <italic>σ</italic>, <italic>n</italic> is the number of states, and <italic>p</italic> is the total number of system parameters. Matrices <italic>U</italic> and <italic>V</italic> contain the eigenvectors of <italic>FF<sup>T</sup></italic> and <italic>F<sup>T</sup>F</italic>, respectively. <italic>U</italic>, Σ, and <italic>V</italic> are ordered according to the magnitude of the singular values. Thus, the first column vector of <italic>U</italic> (and <italic>V</italic>) represents the output (and input) direction with largest amplification. The next most important direction is associated with the second column vector, and so forth. We determine the top three parameters associated with the three greatest input directions in <italic>ν</italic>
<sub>1</sub> and <italic>ν</italic>
<sub>2</sub> as ideal inputs for studying the multiple control input strategy.</p>
<title>Results</title> <p>A 10-state, 38-parameter <italic>Drosophila melanogaster</italic> (fruit fly) circadian model serves as the example system. This stable nonlinear limit cycle oscillator consists of two coupled negative feedback loops that characterize the transcriptional regulation of <italic>period</italic> and <italic>timeless</italic> mRNA and protein dynamics ##REF##9486845##[13]##. <italic>per</italic> and <italic>tim</italic> genes are transcribed in the nucleus, after which their mRNAs are transported into the cytosol where they serve as a template for protein synthesis. The doubly phosphorylated proteins form a heterodimer, PER-TIM, that enters the nucleus and inhibits gene expression, closing the feedback loop. Researchers find that environmental light increases the rate of TIM protein degradation: in this model, light targets the system by magnifying <italic>ν</italic>
<italic><sub>dT</sub></italic>, the doubly phosphorylated TIM protein degradation rate ##REF##9486845##[13]##.</p> <p>The phase response of this model as a function of light is shown via the dash-dotted line in ##FIG##0##Figure 1##. This curve maps the <italic>circadian time</italic> of the entraining stimulus (light pulses) against the resulting change in phase of an organism kept in a free-running environment. The circadian time index repeats every 24 hours with CT0 defining the commencement of dawn and CT12 that of dusk. It is important to note that the magnitude of light-induced phase changes (the quantitative dynamics of the PRC) may vary with respect to the intensity of light. While this model does not account for the complexity of the real network that, for instance, includes additional positive feedback loops ##REF##12581523##[14]##,##REF##11517254##[15]##, it has been experimentally validated ##REF##11517254##[15]## and is widely employed as a reference model ##REF##15340155##[3]##,##REF##11792856##[16]##.</p> <title>Open-Loop Phase Recovery</title> <p>Due to the inherent nonlinear phase response of circadian rhythms when subject to environmental/parametric perturbations, phase recovery dynamics are characterized as a function of the <italic>initial condition</italic> (IC, the circadian time at which control or entrainment begins), and <italic>initial phase difference</italic> (IP, the amount of circadian time to be recovered). To establish a phase resetting set point or upper bound (the maximum amount of time required to recover a given phase difference), we evaluate the open-loop control algorithm, where environmental light∶dark cycles serve as the only mechanism for phase re-entrainment. The phase recovery surface (##FIG##1##Figure 2##) displays the time required for the open-loop case to recover from any possible initial condition and initial phase difference. The asymmetry of the surface may be attributed to the nonlinear effects of light, as characterized by the PRC. The input (light) to output (induced phase shift) mapping of the PRC is seldom symmetric. In <italic>Drosophila melanogaster</italic>, a 15 minute pulse of light has shown to induce up to 3.6 hours of phase advance and 4.2 hours of phase delay ##REF##9486845##[13]##. Recent studies suggest that the change in phase is less sensitive to the duration of the light, and more sensitive to its time-profile ##REF##17502598##[17]##. Phase recovery times (for both open and closed-loop simulations) are evaluated with respect to initial conditions and phase differences discretized at 3 hour intervals. Thus, given the integers <italic>i</italic>,<italic>j</italic> ∈ [0,7], IC = 3<italic>i</italic> and IP = 3<italic>j</italic>.</p> <p>The open-loop entrainment strategy requires at most 183 hours to reset the observed states of the controlled system (cumulative protein complex concentrations) to within 15% of the reference trajectories. Mandating the convergence of state trajectories is a tighter constraint than mandating only phase trajectories, since it incorporates amplitude characteristics. The algorithm, however, may be tuned to consider only strict phase measures. The maximum open-loop recovery time refers to a 9 hour initial phase difference whose control action begins at an initial condition of 15 hours. The initial condition, or start of entrainment, is described with respect to circadian time (CT). Interestingly, there is a stark difference between resetting a 3 to 6 hour initial phase difference versus an 18 to 21 hour initial phase difference (a −6 to −3 hour phase difference). In the former, phase recovers in over 100 hours; in the latter, phase recovers in fewer than 60 hours. Additionally, the open-loop algorithm recovers 9 hour phase differences in a fraction of the time required to correct for smaller phase difference. These properties may be attributed to the nature of the phase response curve and are discussed further in the <xref ref-type="sec" rid="s3">
<italic>Discussion</italic>
</xref> section. Experimental studies in mammalian SCN cells support this asymmetry: Reddy <italic>et al.</italic> show that circadian clock resetting from a 6 hour phase advance (IP6) is accompanied by dissociation of cellular gene expression and may take up to 1 week to recover ##REF##12196553##[18]##. Conversely, resetting a 6 hour phase delay (IP18) is accompanied by coordinated gene expression and requires only 2 days of recovery ##REF##12196553##[18]##. Our simulations support these experimental conclusions as the cumulative protein concentrations in the former case diverge and require several days to converge to the nominal trajectory. In the latter, cumulative protein concentrations oscillate with smaller amplitude until they converge to the nominal trajectory within a couple days. An example of the corresponding simulations is presented in ##SUPPL##0##Figure S1##.</p> <title>Closed-Loop Phase Recovery</title> <p>The MPC algorithm (described in the <italic>Model Predictive Control</italic> section) minimizes the normalized difference between the cumulative protein complex concentration over a prediction horizon of 48 hours, by admitting control action during the first 8 hours of the simulated trajectory. This control action is multiplicative, allowing the algorithm to increase/decrease the nominal parameter by a factor of 2. The control profile defined within the move horizon is updated every 2 hours. Through use of MPC, the re-synchronization rate of the controlled system is increased nearly 3-fold through the control of light, or <italic>ν</italic>
<italic><sub>dT</sub></italic>. Although light serves as a powerful control input, we show that the manipulation of parameters such as transcription and mRNA degradation rates (<italic>ν</italic>
<italic><sub>s</sub></italic> and <italic>ν</italic>
<italic><sub>m</sub></italic>, respectively) may provide more immediate phase resetting. Since we make use of the symmetric version of the mathematical model ##REF##9486845##[13]##, we do not differentiate between <italic>per</italic> or <italic>tim</italic> specific functions. Instead, we assume that the isolated control of <italic>ν</italic>
<italic><sub>sP</sub></italic> is equivalent to the isolated control of <italic>ν</italic>
<italic><sub>sT</sub></italic>, for instance.</p> <p>Parametric sensitivity analysis quantifies the relative change of system behavior with respect to an isolated parametric perturbation. A large sensitivity to a parameter, for instance, suggests that the system's performance is subject to greater change with small variations in the given parameter. We make use of the Fisher Information Matrix (FIM) to evaluate the effect of parametric perturbations on the circadian system's state trajectories ##REF##15695639##[19]##. Investigation of the diagonal values, off diagonal values, and singular value decomposition of the FIM points out the relative order, or rank, of parametric sensitivity measures. This relative ordering highlights sets of control inputs whose manipulation may further reduce phase recovery times. The three greatest diagonal values, for instance, identify the most prominent individual control targets (ranked from most to least sensitive);</p> <p>
<italic>ν</italic>
<italic><sub>s</sub></italic> (the mRNA transcription rate),</p> <p>
<italic>ν</italic>
<italic><sub>m</sub></italic> (the mRNA degradation rate),</p> <p>
<italic>k<sub>s</sub></italic> (the protein translation rate), and</p> <p>
<italic>ν</italic>
<italic><sub>d</sub></italic> (the doubly phosphorylated protein degradation rate).</p> <p>Recall that <italic>ν</italic>
<italic><sub>dT</sub></italic> is the target parameter of environmental light in <italic>Drosophila</italic>. Interestingly, the rate of mRNA transcription is the target of environmental light in <italic>Mus</italic> (via <italic>per</italic> genes) ##REF##12775757##[20]##,##REF##14657377##[21]## and <italic>Neurospora</italic> (via <italic>frq</italic> genes) ##REF##16374510##[22]##. Furthermore, in our previous studies of <italic>Mus</italic> and <italic>Drosophila</italic> circadian networks, mRNA transcription rates were among the most sensitive parameters with respect to both the state- and phase-based sensitivity analysis of two independent network representations ##REF##17158515##[2]##.</p> <p>The greatest off diagonal values identify the most prominent pairs of control targets (ranked accordingly);</p> <p>
<italic>ν</italic>
<italic><sub>s</sub></italic> and <italic>ν</italic>
<italic><sub>m</sub></italic>, and</p> <p>
<italic>ν</italic>
<italic><sub>s</sub></italic> and <italic>K<sub>I</sub></italic> (the threshold constants for repression).</p> <p>Since the manipulation of more than 1 parameter voids the symmetry argument, we target <italic>tim</italic> specific parameters in the implementation of multiple control targets.</p> <p>The greatest input directions of the singular value decomposition identify the most prominent set of three control targets (ranked accordingly);</p> <p>
<italic>ν</italic>
<italic><sub>s</sub></italic>, <italic>ν</italic>
<italic><sub>m</sub></italic>, and <italic>K<sub>I</sub></italic>, and</p> <p>
<italic>ν</italic>
<italic><sub>m</sub></italic>, <italic>ν</italic>
<italic><sub>d</sub></italic>, and <italic>k<sub>2</sub></italic> (the nucleus to cytoplasm rate of transport).</p> <title>Single Target Phase Resetting</title> <p>We investigate the phase recovery dynamics corresponding to four independent isolated control inputs with respect to the initial condition and initial phase difference (##FIG##2##Figure 3##). Results show that control targets identified via sensitivity analysis (##FIG##2##Figure 3(A)–3(C)##) serve as more effective re-entrainment factors than light (##FIG##2##Figure 3(D)##). More specifically, the maximum recovery time corresponding to a control input of <italic>ν<sub>s</sub></italic> is 44 hours (at IC9 and IP12/IP15), <italic>ν</italic>
<italic><sub>m</sub></italic> is 50 hours (at IC21 and IP15), <italic>k<sub>s</sub></italic> is 59 hours (at IC12 and IP15), and <italic>ν</italic>
<italic><sub>d</sub></italic> (the light target) is 60 hours (at IC12 and IP15, or IC9 and IP12). The control profiles and state response dynamics relating to the phase recovery of IC9 and IP12 are provided in ##SUPPL##1##Figure S2## and ##SUPPL##2##Figure S3##.</p> <p>There is a subtle similarity among the single-input phase recovery data; namely, the sudden drop in recovery time with respect to the initial condition for initial phase differences of 0 to 15 hours. We attribute this steep recovery gradient to the PRC as it depicts a greater region of phase delay than it does a phase advance. For this reason, it is more beneficial if the organism delays its phase to recover from a 12 hour initial phase difference. Furthermore, recall that a phase delay is incurred if the organism is to receive a photic input in the late evening hours. Hence, recovering from a phase difference via a set of delaying control inputs is most efficient if control action begins around the late subjective evening. Thus, if we observe phase resetting behavior corresponding to a small phase difference (such that the subjective day of the controlled system and reference are similar), we expect it to have the shortest recovery time near an initial condition of 12 hours, or dusk (##FIG##2##Figure 3(D)##). Interestingly, each of the control inputs exhibits this property. We attribute this similarity to the unique PRC of each control input (##FIG##0##Figure 1##).</p> <title>Dual target phase resetting</title> <p>Allowing the MPC algorithm to manipulate two variables simultaneously provides more immediate phase resetting since the controller has greater flexibility. For instance, the simultaneous use of <italic>ν<sub>s</sub></italic> and <italic>ν</italic>
<italic><sub>m</sub></italic> requires a maximum phase recovery of 43 hours to recover a 12 hour initial phase difference entrained from an initial condition of 6 or 9 hours (##FIG##3##Figure 4(A)##). Similarly, the simultaneous control of <italic>ν</italic>
<italic><sub>s</sub></italic> and <italic>K<sub>I</sub></italic> requires a maximum phase recovery of 46 hours to recover a 15 hour initial phase difference (##FIG##3##Figure 4B##).</p> <p>Given that the common input is <italic>ν</italic>
<italic><sub>s</sub></italic>, we expect the dual control input phase recovery dynamics to be just as good (if not better) than the results generated from the single <italic>ν</italic>
<italic><sub>s</sub></italic> input algorithm. Although the dual control input strategy provides similar maximum phase recovery times, the greatest recovery time “plateau” is smaller. Therefore, the dual <italic>ν</italic>
<italic><sub>s</sub></italic> and <italic>ν<sub>m</sub></italic> input strategy is more effective at recovering an initial phase differences of 15 hours from IC6, while the <italic>ν</italic>
<italic><sub>s</sub></italic> and <italic>K<sub>I</sub></italic> pair is more effective at recovering a 12 hour initial phase difference from IC9 (compare ##FIG##2##Figure 3(A)## to ##FIG##3##Figure 4(A)–4(B)##). We would even argue that the recovery times associated with the dual input strategy may lessen if the genetic algorithm based optimizer were run over a greater number of generations. We limit the number of generations in the genetic algorithm – 15 for the single input, 75 for the dual input, and 250 for the triple input – to reflect the limited resources and time constraints evident in real world applications.</p> <title>Triple target phase resetting</title> <p>Just as the dual input case, we expect the triple input strategy to recover phase just as effectively as the single input strategy. The <italic>ν<sub>s</sub></italic>, <italic>ν</italic>
<italic><sub>m</sub></italic>, and <italic>K<sub>I</sub></italic> input strategy requires at most 39 hours to recover an initial phase difference of 12 hours at IC9 (##FIG##3##Figure 4(C)##). Interestingly, the maximum recovery time corresponding to the use of <italic>ν</italic>
<italic><sub>m</sub></italic>, <italic>ν</italic>
<italic><sub>d</sub></italic>, and <italic>k<sub>2</sub></italic> as simultaneous control inputs is 59 hours to recover a 3 hour phase difference at IC12 (##FIG##3##Figure 4(D)##). We attribute this abnormally high recovery time to the possibility of numerical errors associated with the optimization algorithm since <italic>ν</italic>
<italic><sub>m</sub></italic> requires no more than 30 hours to recover from a 3 hour phase difference, while <italic>ν</italic>
<italic><sub>d</sub></italic> requires no more than 54 hours. If we omit this data point as an outlier, this triple input strategy requires 42 hours to recover a 12 hour phase difference beginning at IC6. Assuming abundant computational resources and time, the triple input strategy may further outperform the dual input strategy since the MPC algorithm acquires greater flexibility (a greater number of control options) with each additional target. Each of these control inputs produces a unique PRC that allows the algorithm to further manipulate the set of targets such that the combination may yield a phase delay or advance at any time of the circadian day. In the case of a single light (or <italic>ν</italic>
<italic><sub>d</sub></italic> target) input, for instance, the algorithm must wait for the subjective morning to force a phase advance, or the subjective night to force a delay. The advantage gained through additional control targets, however, is not clear. Given the finite horizon over which the algorithm optimizes phase synchrony, in addition to the nonlinear response of the model, we can not expect a monotonic improvement of phase recovery dynamics with an increase in the number of manipulated variables. For instance, the algorithm may choose a sequence of multiple inputs that yields lower cost in the short term (as compared to a single input) with a greater cost in the long term, leading to a point of no return. This scenario may also attribute to the 59 hour recovery observed in ##FIG##3##Figure 4(D)##.</p> <title>Short and Long Period Mutants</title> <p>Mutant phenotypes of the circadian oscillator represent cases in which nominal light∶dark cycles are unable to maintain synchrony. For this reason, the MPC tuning parameters must be re-evaluated according to this phase resetting problem. In wild-type, for instance, we can afford to be more aggressive with control penalties since nominal light∶dark cycles (or, no control) will eventually entrain the system. In mutants, the weights used to penalize the state error and control inputs prove to be more influential since nominal light∶dark cycles will not entrain the system. Therefore, we set both the move and prediction horizon to 24 hours and reduce the penalty of state error and control to ones. To counter the computational expense incurred with a longer move horizon, we set the time step to 4 hours. Through MPC, we identify a more suitable light∶dark cycle that synchronizes organisms exhibiting abnormally short and long free-running periods (22 and 27 hours, respectively, as shown in ##FIG##4##Figure 5##). Determining the complete range of entrainment (which is likely wider than the 22 to 27 hour period) is non-trivial. In a previous study, we found that (i) the predicted range of entrainment may be very sensitive to the employed performance metric, and (ii) the control/light input strength may also play a dominant role in defining the bounds of this range ##UREF##6##[23]##.</p> <p>Given that the PRC characterizing the behavior of <italic>Drosophila melanogaster</italic> consists of phase delays during the late subjective evening, we expect short-period mutant phenotypes to require bright light after subjective dusk. Similarly, we expect long-period mutant phenotypes to require bright light in the early subjective morning to advance the cycle. Our results confirm this hypothesis. In ##FIG##1##Figure 2##, we demonstrate how bright light, admitted during the environmental night, resets the phase of short-period mutants such that it matches that of its environment. Given that the controlled system is 2 hours short, the occurrence of light during the night overlaps with the advance region of the system's PRC. Similarly, the onset of bright light at dawn overlaps with the delay region of long-period mutant PRCs (##FIG##1##Figure 2##). Our ability to maintain appropriate phase relationships between mutant phenotypes (models characterized by non-nominal parameters) and the environment (the nominal case) further proves the robustness of the algorithm despite model mismatch.</p>
<title>Discussion</title> <title>Circadian Phase Response</title> <p>As implied by the PRC (##FIG##0##Figure 1##), a 3 hour phase difference may be recovered immediately through admission of light at CT15. Hence, for open-loop control action to be most effective, environmental daylight should occur during the controlled system's subjective night (at CT15). In cases with small initial phase difference (such that the subject's internal time is nearly equal to environmental time), however, daylight begins entrainment once the subjective day is around CT12, by inducing small phase delays. This delay reduces the overlap between environmental daylight and the subjective night since re-entrainment of the initial phase difference began before subjective night. The opposite occurs with small <italic>negative</italic> phase differences, where an 18 hour (or −6 hour) phase difference may be recovered via a light pulse admitted at CT21. In this case, environmental daylight affects the controlled system at the start of day while it has not yet begun entrainment, maximizing the phase advancing effect of light. For this reason, open-loop entrainment via phase advances requires less recovery time despite the fact that a single pulse of light may induce a greater phase delay than advance.</p> <p>More generally, we find that any given initial phase difference is more readily recovered if open-loop entrainment begins between CT0 and CT9; the rate of re-entrainment depends on the initial condition. To correct initial phase differences of 0 to 9 hours (by inducing a phase delay), daylight is most effective at the end of the day, suggesting greater performance if the algorithm were to begin control action around CT6. To correct initial phase differences of 0 to −6 hours (by inducing phase advances), daylight is most effective at the start of the day, suggesting greater performance if the algorithm were to begin around CT0. In the former case, daylight overlaps with the delay region of the subject's PRC, while in the latter it overlaps with the advance region. Resetting an initial condition of 12 to 15 hours, however, presents an interesting control dilemma as environmental daylight may induce both a phase delay and phase advance. For this reason, the open-loop control algorithm requires several days to correct for such phase differences. If light were accessible to entrain the system continuously throughout the day and night (in other words, if we were to <italic>close the loop</italic>), phase recovery dynamics would be less extreme since phase resetting would rely less on the initial condition.</p> <p>Additional phase resetting properties may be inferred through investigation of the simulated PRCs. For instance, in the single input case, <italic>ν</italic>
<italic><sub>s</sub></italic> and <italic>k<sub>s</sub></italic> exhibit similar recovery dynamics with the exception that <italic>ν</italic>
<italic><sub>s</sub></italic> is more effective at resetting initial phase differences of 15 to 21 hours. This quality may be associated with the fact that manipulating <italic>k<sub>s</sub></italic> exhibits a strikingly similar phase response as <italic>ν</italic>
<italic><sub>s</sub></italic> where their input to output mapping is shifted by about 5 hours (##FIG##0##Figure 1##). This similarity may be attributed to the fact <italic>k<sub>s</sub></italic> and <italic>ν</italic>
<italic><sub>s</sub></italic> are directly involved with the irreversible production, and transcriptional/translational regulation, of clock-specific genes/proteins. Additionally, the “active” region of the <italic>ν</italic>
<italic><sub>s</sub></italic> and <italic>ν<sub>m</sub></italic> PRCs are wider than those of <italic>k<sub>s</sub></italic> and <italic>ν</italic>
<italic><sub>d</sub></italic> (or, their dead zones are shorter than those of <italic>k<sub>s</sub></italic> and <italic>ν</italic>
<italic><sub>d</sub></italic>), suggesting that their perturbation-induced phase shifts are accessible throughout a greater portion of the circadian day.</p> <title>Minimizing Control Action</title> <p>Of the single control input results, the manipulation of <italic>ν</italic>
<italic><sub>s</sub></italic>, identified as the most sensitive parameter, provides the shortest phase recovery times. Despite these results, <italic>ν</italic>
<italic><sub>d</sub></italic> or light-based control is most efficient. In ##FIG##5##Figure 6##, we relate the cumulative control input (a unitless measure that integrates the multiplicative control target action) to the convergence of phase via the PER-TIM complex state error. The data shown reflects the recovery of an initial phase difference of 15 hours from IC12. Analyzing this relationship may provide a basis from which the pharmaceutical industry might select one drug over another. If two different drug targets demonstrate similar response, the one that requires the least number of doses should be admitted, minimizing cost and the potential for drug related side-effects. Moreover, if the symptoms of illness are more severe than the potential for side effect, the drug that minimizes the state error may be preferred over others. The assessment of system convergence and the corresponding admitted control is key to the identification and application of control targets.</p> <title>Circadian Alignment and Illness</title> <p>In our modern “24/7” work world, social and commercial pressures often oppose our natural circadian timekeeping, causing a source of circadian stress that may lead to chronic illnesses such as cardiovascular disease and cancer ##UREF##7##[24]##. Numerous studies seem to show the effect of circadian rhythms on processes such as cell proliferation and apoptosis that eventually lead to proper growth control ##REF##14633665##[25]##–##REF##16809488##[27]##. For instance, components of the cell cycle that dictate the G1-S and G2-M transition phase have been associated with circadian transcriptional regulation ##REF##14523398##[28]##,##REF##12934012##[29]##. Also in certain conditions, cancer can be a direct consequence of the absence of the circadian regulation ##REF##14633665##[25]##,##REF##12724733##[26]##,##REF##12491517##[30]##. A review of circadian related clinical disorders describes how mutations in some clock genes are associated with alcoholism, sleeping disorders, hypertension, and morbidity ##UREF##7##[24]##,##REF##16077156##[31]##. Most commonly, poor circadian regulation leads to advanced sleep phase syndrome, delayed sleep phase syndrome, non-24-hour sleep-wake syndrome, and irregular sleep-wake pattern ##UREF##8##[32]##. In each of these cases, poor circadian phase resetting may be achieved through the systematic admission of controlled light pulses.</p> <p>Assuming we have access to drugs that specifically target circadian genes, we can identify the targets whose manipulation yields the most effective and immediate response through investigation of each control's phase dynamics (as shown in ##FIG##0##Figure 1##). Or, it is possible to minimize the use of control and choose targets that require the least number of doses. We may also tailor the MPC algorithm to correct phase more readily through simultaneous manipulation of multiple control targets. Even further, we may reduce the computational expense by enumerating the control solutions over a grid in the solution space (light magnitude as a function of time), and choosing the optimal control sequence via an exhaustive search. The algorithm approaches a globally optimal solution as the total possible quantization steps of the control input increases. We tested the efficacy of the algorithm with respect to a quantization of 2, 4, 8, and 16 steps ##UREF##5##[12]##. Results suggest that the shorter recovery time associated with the finer-grid enumeration may not outweigh the increase in computation time. Therefore, we may dramatically reduce computational expense by investigating control solutions for as few as 2 possible control values.</p> <p>Our methods show great promise for use in the pharmaceutical industry as our theoretical phase entrainment of mutant phenotypes demonstrates the robustness of the algorithm in the presence of model mismatch. This robustness alleviates concerns in the pharmaceutical industry to tailor mathematical representations of bio-chemical pathways to individual people.</p> <title>Mammalian Studies</title> <p>The study of controlled light pulses as a means of correcting phase is a common area of interest. Studies have shown that humans are much more sensitive to light than initially suspected since room light can significantly reset the phase of the human circadian clock ##UREF##9##[33]##,##REF##8596632##[34]##. Furthermore, the admission of morning light has been considered as an antidepressant by realigning the internal clock with the environment ##REF##16648247##[35]##.</p> <p>Additional studies suggest that the human circadian clock mechanism functions similarly to those of other mammals ##REF##8596632##[34]##. This similarity may be attributed to shape/amplitude characteristics of their respective phase response curves. Humans show phase-delay shifts of up to 3.6 hours and phase-advance shifts of up to 2.01 hours (with respect to a 6.7 hour pulse of bright light) ##REF##12717008##[36]##, which is both quantitatively and qualitatively similar to other mammalian species. This parallel motivates the experimental application of controlled light pulses for phase resetting in mammals. We have taken this first step by assessing the efficacy and computational utility of model predictive control as applied to a detailed 71-state <italic>Mus musculus</italic> circadian model ##UREF##10##[37]##. Furthermore, melatonin has proven to be a key circadian phase resetting agent for totally blind people who cannot synchronize to environmental day∶night cycles (or do so at an abnormal time) ##REF##16648247##[35]##. Therefore, melatonin may be used individually (in cases to treat the totally blind), or in combination with light to provide more effective phase resetting.</p> <p>Therapies designed to alleviate circadian load would have an important impact on morbidity and mortality across the developed world. Aside from correcting mutant phenotypes, phase resetting would increase performance in many healthy, or wild-type, cases such as frequent flyers avoiding jet-lag or astronauts maintaining a rigorous schedule during space exploration ##REF##17502598##[17]##. The real-time application of the proposed algorithm, however, may be a major issue; in practice, it will not be feasible to collect the corresponding protein concentration data at the molecular level. However, behavioral and/or physiological parameters that are controlled by (and correlated with) the circadian clock's dynamics are easily accessible. Such data may include actograms such as wheel running data for rodents ##UREF##3##[6]##. Hence, a missing link in the current work concerns the development of corresponding (non-linear) state estimators for reconstructing the molecular dynamics. Given the discrete nature of MPC (sampling every 4 hours), the proposed strategy is feasible in practice since sampling rates of such physiological circadian markers may be much higher.</p>
<p>Conceived and designed the experiments: NB JS FD. Performed the experiments: NB. Analyzed the data: NB JS FD. Wrote the paper: NB JS FD. Principal Investigator: FD.</p> <p>Circadian entrainment is necessary for rhythmic physiological functions to be appropriately timed over the 24-hour day. Disruption of circadian rhythms has been associated with sleep and neuro-behavioral impairments as well as cancer. To date, light is widely accepted to be the most powerful circadian synchronizer, motivating its use as a key control input for phase resetting. Through sensitivity analysis, we identify additional control targets whose individual and simultaneous manipulation (via a model predictive control algorithm) out-perform the open-loop light-based phase recovery dynamics by nearly 3-fold. We further demonstrate the robustness of phase resetting by synchronizing short- and long-period mutant phenotypes to the 24-hour environment; the control algorithm is robust in the presence of model mismatch. These studies prove the efficacy and immediate application of model predictive control in experimental studies and medicine. In particular, maintaining proper circadian regulation may significantly decrease the chance of acquiring chronic illness.</p> <title>Author Summary</title> <p>The robust timing, or phase, of the circadian clock is critical in directing and synchronizing molecular, cellular, and organismal behaviors. The clock's failure to maintain precision and adaption is associated with sleeping disorders, depression, and cancer. To better study and control the timing of circadian rhythms, we make use of systems theoretic tools such as sensitivity analysis and model predictive control (MPC). Sensitivity analysis is used to identify key driving mechanisms without having to fully understand or investigate the detailed mechanistic interconnections of the large complex circadian network. Contrary to intuition, sensitivity analysis of the circadian model highlights several <italic>non</italic>-photic control inputs (such as transcriptional regulation) that outperform light-based circadian phase resetting – light is known to accelerate protein degradation. Aside from targeting individual parameters as control inputs, our methods identify combinations of control targets that may further the efficiency of entrainment. We compare the phase resetting performance of our MPC algorithm among cases involving individual and multiple simultaneous control targets (in wild-type simulations). We then tailor the algorithm to correct continuously the phase mismatch that occurs in short and long period mutant phenotypes. Through use of the presented tools, our algorithm is robust in the presence of model mismatch and outperforms the natural <italic>in silico</italic> sun-cycle–based phase recovery strategy by nearly 3-fold.</p>
<title>Supporting Information</title>
<fig id="pcbi-1000104-g001" position="float"><object-id pub-id-type="doi">10.1371/journal.pcbi.1000104.g001</object-id><label>Figure 1</label><caption><title>Circadian phase response behavior.</title><p>Phase response curves traditionally characterize the light pulse to induced phase mapping of the input admitted to a free-running circadian oscillator. Here, phase response dynamics of the four system parameters exhibiting greatest state sensitivity is depicted: <italic>ν<sub>s</sub></italic> (mRNA transcription), <italic>ν<sub>m</sub></italic> (mRNA degradation), <italic>k<sub>s</sub></italic> (protein translation), and <italic>ν<sub>d</sub></italic> (protein degradation). The x-axis denotes the time at which the 2 hour pulse is given (where CT0 reflects dawn and CT12 dusk), and the y-axis describes the induced phase shift. A positive shift reflects a phase advance. Since light targets TIM specific protein degradation, <italic>ν<sub>dT</sub></italic>, the light-based PRC of the <italic>Drosophila</italic> model is represented via the dash-dotted line.</p></caption></fig> <fig id="pcbi-1000104-g002" position="float"><object-id pub-id-type="doi">10.1371/journal.pcbi.1000104.g002</object-id><label>Figure 2</label><caption><title>Natural phase entrainment.</title><p>Open-loop phase resetting dynamics are plotted as a function of the initial phase difference (x-axis) and the initial condition (y-axis). The intensity of the color reflects the amount of time required to reset a given phase via the light∶dark cycles calibrated to the initial condition: the lighter the color, the longer the recovery time. The mapping of color intensity to phase recovery times (in hours) is shown in the vertical color bar.</p></caption></fig> <fig id="pcbi-1000104-g003" position="float"><object-id pub-id-type="doi">10.1371/journal.pcbi.1000104.g003</object-id><label>Figure 3</label><caption><title>Single input control.</title><p>Closed-loop phase resetting dynamics for single control targets (ordered according to their relative sensitivity) are described as a function of the initial phase difference (x-axis) and the initial condition (y-axis). The intensity of the color gradient reflects the amount of time required to recover from the given control conditions: the lighter the color, the longer the phase recovery. Each color bar is calibrated according to a minimum recovery time of 0 hours and maximum of 60 hours. (A) <italic>ν</italic>
<italic><sub>s</sub></italic> Single Control Target (B) <italic>ν</italic>
<italic><sub>m</sub></italic> Single Control Target (C) <italic>k<sub>s</sub></italic> Single Control Target (D) <italic>ν</italic>
<italic><sub>d</sub></italic> Single Control Target.</p></caption></fig> <fig id="pcbi-1000104-g004" position="float"><object-id pub-id-type="doi">10.1371/journal.pcbi.1000104.g004</object-id><label>Figure 4</label><caption><title>Multiple input control.</title><p>Closed-loop phase resetting dynamics for dual ((A) and (B)), and triple ((C) and (D)) control targets are shown with respect to the initial phase difference (x-axis) and the initial condition (y-axis). The phase recovery time is denoted by the intensity of the color at each given data point: the lighter the color, the longer the recovery time. The mapping of color intensity to the recovery time (in hours) is reflected in the color bar. Each color bar is calibrated according to a minimum recovery time of 0 hours and maximum of 60 hours. (A) <italic>ν</italic>
<italic><sub>s</sub></italic> and <italic>ν</italic>
<italic><sub>m</sub></italic> Dual Control Targets (B) <italic>ν</italic>
<italic><sub>s</sub></italic> and <italic>K<sub>I</sub></italic> Dual Control Targets (C) <italic>ν</italic>
<italic><sub>s</sub></italic>, <italic>ν</italic>
<italic><sub>m</sub></italic>, and <italic>K<sub>I</sub></italic> Triple Control Targets (D) <italic>ν</italic>
<italic><sub>m</sub></italic>, <italic>ν</italic>
<italic><sub>d</sub></italic>, and <italic>k</italic>
<sub>2</sub> Triple Control Targets.</p></caption></fig> <fig id="pcbi-1000104-g005" position="float"><object-id pub-id-type="doi">10.1371/journal.pcbi.1000104.g005</object-id><label>Figure 5</label><caption><title>Model mis-match.</title><p>Continuous phase resetting for the short-period mutant phenotype (dotted lines) and long-period mutant phenotype (dashed lines) are depicted with respect to an initial condition of 0 and 12 hours. Upper subplots describe the observed state trajectory (cumulative PER-TIM protein complex concentrations) as a function of controlled light pulses, shown in the lower subplots. The nominal response (denoted by solid lines) is entrained via regular 24 hour light∶dark cycles. As expected, short-period mutants reset via daily light pulses that occur during the subjective night, forcing daily phase delays. Long period mutants reset via daily light pulses that occur during the subjective morning, forcing phase advances. (A) Mutant Response at IC0 (B) Mutant Response at IC12.</p></caption></fig> <fig id="pcbi-1000104-g006" position="float"><object-id pub-id-type="doi">10.1371/journal.pcbi.1000104.g006</object-id><label>Figure 6</label><caption><title>The cost of control.</title><p>The convergence of state dynamics (from an initial phase difference of 15 hours at IC12) is plotted (top) as a function time. The corresponding cumulative control input (determined by integrating the value of each multiplicative control target over time) is described in the lower subplot. The light-induced target, <italic>ν</italic>
<italic><sub>d</sub></italic>, shown via the dash-dotted line, requires the greatest amount of multiplicative control while exhibiting the greatest amount of state error. Conversely, the <italic>k<sub>s</sub></italic> target corresponds to the least amount of state error and requires less admitted control.</p></caption></fig>
<disp-formula></disp-formula> <inline-formula></inline-formula> <inline-formula></inline-formula> <inline-formula></inline-formula> <inline-formula></inline-formula> <inline-formula></inline-formula> <inline-formula></inline-formula> <disp-formula></disp-formula> <inline-formula></inline-formula> <inline-formula></inline-formula> <disp-formula></disp-formula> <inline-formula></inline-formula> <inline-formula></inline-formula>
<supplementary-material content-type="local-data" id="pcbi.1000104.s001"><label>Figure S1</label><caption><p>Open loop phase resetting response at IC9. Phase resetting dynamics for an initial phase difference of 18 hours (or −6 hours) is shown in the blue dotted trajectory; those pertaining to IP6 are reflected in the red dashed line. The nominal protein concentration dynamics are depicted in the solid black line, while environmental sun cycles are shown in the black dotted square wave. The magnitude of the square wave oscillates between 1 and 2 and does not correspond to the y-axis of the figure.</p><p>(0.02 MB EPS)</p></caption></supplementary-material> <supplementary-material content-type="local-data" id="pcbi.1000104.s002"><label>Figure S2</label><caption><p>Single control input to output response for IC9. Phase resetting dynamics (upper subplot) are depicted as a function of the individual control input profiles (lower subplot). Nominal, or pre-entrained, circadian dynamics are shown in solid black. The blue dotted lines reflect phase resetting with respect to <italic>ν</italic>s, while the red dashed lines reflect those of <italic>ν</italic>m. Although phase resetting, or state convergence, among the four different control variable occurs at similar hours, both the state dynamics and control profiles for each variable are significantly different.</p><p>(0.04 MB EPS)</p></caption></supplementary-material> <supplementary-material content-type="local-data" id="pcbi.1000104.s003"><label>Figure S3</label><caption><p>Single control input to output response for IP12. Phase resetting dynamics (upper subplot) are depicted as a function of the individual control input profiles (lower subplot). Nominal, or pre-entrained, circadian dynamics are shown in solid black. The blue dotted lines reflect phase resetting with respect to ks, while the red dashed lines reflect those of <italic>ν</italic>d. Although phase resetting, or state convergence, among the four different control variable occurs at similar hours, both the state dynamics and control profiles for each variable are significantly different.</p><p>(0.03 MB EPS)</p></caption></supplementary-material>
<fn-group><fn fn-type="COI-statement"><p>The authors have declared that no competing interests exist.</p></fn><fn fn-type="financial-disclosure"><p>This project was supported in part by the ICB, DAAD19-03-D-0004; NIH, GM078993; The Research Participation Program between the US DOE and AFRL/HEP.</p></fn></fn-group>
<graphic xlink:href="pcbi.1000104.g001"/> <graphic xlink:href="pcbi.1000104.g002"/> <graphic xlink:href="pcbi.1000104.g003"/> <graphic xlink:href="pcbi.1000104.g004"/> <graphic xlink:href="pcbi.1000104.g005"/> <graphic xlink:href="pcbi.1000104.g006"/> <graphic xlink:href="pcbi.1000104.e001.jpg" mimetype="image" position="float"/> <inline-graphic xlink:href="pcbi.1000104.e002.jpg" mimetype="image"/> <inline-graphic xlink:href="pcbi.1000104.e003.jpg" mimetype="image"/> <inline-graphic xlink:href="pcbi.1000104.e004.jpg" mimetype="image"/> <inline-graphic xlink:href="pcbi.1000104.e005.jpg" mimetype="image"/> <inline-graphic xlink:href="pcbi.1000104.e006.jpg" mimetype="image"/> <inline-graphic xlink:href="pcbi.1000104.e007.jpg" mimetype="image"/> <graphic xlink:href="pcbi.1000104.e008.jpg" mimetype="image" position="float"/> <inline-graphic xlink:href="pcbi.1000104.e009.jpg" mimetype="image"/> <inline-graphic xlink:href="pcbi.1000104.e010.jpg" mimetype="image"/> <graphic xlink:href="pcbi.1000104.e011.jpg" mimetype="image" position="float"/> <inline-graphic xlink:href="pcbi.1000104.e012.jpg" mimetype="image"/> <inline-graphic xlink:href="pcbi.1000104.e013.jpg" mimetype="image"/>
<media xlink:href="pcbi.1000104.s001.eps"><caption><p>Click here for additional data file.</p></caption></media> <media xlink:href="pcbi.1000104.s002.eps"><caption><p>Click here for additional data file.</p></caption></media> <media xlink:href="pcbi.1000104.s003.eps"><caption><p>Click here for additional data file.</p></caption></media>
[{"label": ["1"], "element-citation": ["\n"], "surname": ["Zeilinger", "Farre", "Taylor", "Kay", "Doyle"], "given-names": ["M", "E", "S", "S", "F"], "suffix": ["III"], "year": ["2006"], "article-title": ["A novel computational model of the circadian clock in Arabidopsis that incorporates prr7 and prr9."], "source": ["Mol Syst Biol"], "volume": ["2"], "fpage": ["1"], "lpage": ["13"]}, {"label": ["4"], "element-citation": ["\n"], "surname": ["Bagheri", "Stelling", "Doyle"], "given-names": ["N", "J", "F"], "suffix": ["III"], "year": ["2007"], "article-title": ["Circadian phase entrainment via nonlinear model predictive control."], "source": ["Int J Robust Nonlinear Control"], "comment": ["doi:10.1002/rnc.1209"]}, {"label": ["5"], "element-citation": ["\n"], "surname": ["Mott", "Mollicone", "van Wollen", "Huzmezan"], "given-names": ["C", "D", "M", "M"], "year": ["2003"], "article-title": ["Modifying the human circadian pacemaker using model based predictive control."], "publisher-loc": ["Denver, CO"], "publisher-name": ["American Control Conference"], "fpage": ["453"], "lpage": ["458"]}, {"label": ["6"], "element-citation": ["\n"], "surname": ["Daan", "Pittendrigh"], "given-names": ["S", "CS"], "year": ["1976"], "article-title": ["A functional analysis of circadian pacemakers in nocturnal rodents. II. The variability of phase response curves."], "source": ["J Comp Physiol"], "volume": ["106"], "fpage": ["253"], "lpage": ["266"]}, {"label": ["10"], "element-citation": ["\n"], "surname": ["Slaby", "Sager", "Shaik", "Kummer", "Lebiedz"], "given-names": ["O", "S", "OS", "U", "D"], "year": ["2007"], "article-title": ["Optimal control of self-organized dynamics in cellular signal transduction."], "source": ["Math Comp Model Dyn"], "volume": ["13"], "fpage": ["487"], "lpage": ["502"]}, {"label": ["12"], "element-citation": ["\n"], "surname": ["Taylor", "Bagheri", "Meeker", "Petzold", "Doyle"], "given-names": ["SR", "N", "K", "LR", "FJ"], "suffix": ["III"], "year": ["2007"], "article-title": ["Robust timekeeping in circadian networks: From genes to cells."], "comment": ["In: Proceedings of the Foundations of Systems Biology in Engineering; 9\u201312 September 2007; Stuttgart, Germany. FOSBE 2007"]}, {"label": ["23"], "element-citation": ["\n"], "surname": ["Doyle", "Stelling"], "given-names": ["FJ", "J"], "suffix": ["III"], "year": ["2005"], "article-title": ["Robust performance in biophysical networks."], "publisher-loc": ["Prague, Czech"], "publisher-name": ["IFAC World Congress"]}, {"label": ["24"], "element-citation": ["\n"], "surname": ["Hastings", "Reddy", "Maywood"], "given-names": ["MH", "AB", "ES"], "year": ["2003"], "article-title": ["A clockwork web: Circadian timing in brain and periphery in health and disease."], "source": ["Nat Neurosci"], "volume": ["4"], "fpage": ["649"], "lpage": ["661"]}, {"label": ["32"], "element-citation": ["\n"], "surname": ["Lamberg"], "given-names": ["L"], "year": ["1994"], "source": ["Bodyrhythms: chronobiology and peak performance"], "publisher-loc": ["New York, NY"], "publisher-name": ["William Morrow and Company"]}, {"label": ["33"], "element-citation": ["\n"], "surname": ["Gronfier", "Wright", "Kraunauer", "Jewett"], "given-names": ["C", "KP", "RE", "ME"], "suffix": ["Jr"], "year": ["2004"], "article-title": ["Efficacy of a single sequence of intermittent bright light pulses for delaying circadian phase in humans."], "source": ["Am J Physiol Endocrinol Metab"], "volume": ["287"], "fpage": ["174"], "lpage": ["181"]}, {"label": ["37"], "element-citation": ["\n"], "surname": ["Bagheri", "Taylor", "Meeker", "Petzold", "Doyle"], "given-names": ["N", "SR", "K", "LR", "FJ"], "suffix": ["III"], "year": ["2008"], "article-title": ["Synchrony and entrainment properties of robust circadian oscillators."], "source": ["J R Soc Interface"], "comment": ["doi:10.1098/rsif.2008.0045.focus"]}, {"label": ["38"], "element-citation": ["\n"], "surname": ["Morari", "Lee"], "given-names": ["M", "JH"], "year": ["1999"], "article-title": ["Model predictive control: Past, present and future."], "source": ["Comp Chem Eng"], "volume": ["23"], "fpage": ["667"], "lpage": ["682"]}, {"label": ["40"], "element-citation": ["\n"], "surname": ["Beyer", "Schwefel"], "given-names": ["HG", "HP"], "year": ["2002"], "article-title": ["Evolution strategies."], "source": ["Nat Comput"], "volume": ["1"], "fpage": ["3"], "lpage": ["52"]}, {"label": ["41"], "element-citation": ["\n"], "surname": ["Goldberg"], "given-names": ["DE"], "year": ["1989"], "source": ["Genetic Algorithms in Search, Optimization, and Machine Learning"], "publisher-loc": ["Reading, MA"], "publisher-name": ["Addison-Wesley Publishing Company"]}, {"label": ["42"], "element-citation": ["\n"], "surname": ["Michalewicz"], "given-names": ["Z"], "year": ["1996"], "source": ["Genetic Algorithms+Data Structures\u200a=\u200aEvolution Programs"], "publisher-loc": ["New York, NY"], "publisher-name": ["Springer-Verlag"], "comment": ["3rd edition"]}, {"label": ["43"], "element-citation": ["\n"], "surname": ["Skogestad", "Postlethwaite"], "given-names": ["S", "I"], "year": ["1996"], "source": ["Multivariable Feedback Control."], "publisher-name": ["John Wiley & Sons"]}]
PLoS Comput Biol. 2008 Jul 4; 4(7):e1000104