Optimal Control of Oscillation Timing and Entrainment Using Large Magnitude Inputs: An Adaptive Phase-Amplitude-Coordinate-Based Approach
OOptimal Control of Oscillation Timing and Entrainment UsingLarge Magnitude Inputs: An AdaptivePhase-Amplitude-Coordinate-Based Approach
Dan Wilson Department of Electrical Engineering and Computer Science, University of Tennessee,Knoxville, TN 37996, USAFebruary 10, 2021
Abstract
Given the high dimensionality and underlying complexity of many oscillatory dynamicalsystems, phase reduction is often an imperative first step in control applications where oscillationtiming and entrainment are of interest. Unfortunately, most phase reduction frameworks placerestrictive limitations on the magnitude of allowable inputs, limiting the practical utility of theresulting phase reduced models in many situations. In this work, motivated by the search forcontrol strategies to hasten recovery from jet-lag caused by rapid travel through multiple timezones, the efficacy of the recently developed adaptive phase-amplitude reduction is considered formanipulating oscillation timing in the presence of a large magnitude entraining stimulus. Theadaptive phase-amplitude reduced equations allow for a numerically tractable optimal controlformulation and the associated optimal stimuli significantly outperform those resulting fromfrom previously proposed optimal control formulations. Additionally, a data-driven techniqueto identify the necessary terms of the adaptive phase-amplitude reduction is proposed andvalidated using a model describing the aggregate oscillations of a large population of coupledlimit cycle oscillators. Such data-driven model reduction algorithms are essential in situationswhere the underlying model equations are either unreliable or unavailable.
Mathematical analysis and control of oscillatory dynamical systems is a widely studied problemwith relevant applications to neurological brain rhythms [37], [20], [26], circadian physiology [11],[3], [48], and various other physical and chemical systems [50], [13], [45], [62], [66]. Given the highdimensionality and sheer complexity of many oscillatory dynamical systems, model reduction isoften a necessary first step in their analysis. While many reduction strategies have been developedfor limit cycle oscillators, phase reduction [64], [14], [27] is one of the most widely applied techniqueswhich allows for oscillatory systems of the form˙ x = F ( x ) + U ( t ) , (1)where x ∈ R N is the system state, F describes the nominal dynamics, and U is an exogenous inputto be analyzed according to ˙ θ = ω + Z ( θ ) · U ( t ) , (2)1 a r X i v : . [ m a t h . D S ] F e b here the phase θ ∈ [0 , π ) characterizes the state of the oscillator in reference to a limit cycle, ω ( p ) = 2 π/T is the natural frequency where T is the nominal, unperturbed period, Z ( θ ) is thephase response curve (i.e., the gradient of the phase coordinate evaluated on the limit cycle),and the ‘dot’ denotes the dot product. Phase reduction provides a tremendous decrease in thedimensionality of a given limit cycle oscillator, allowing for the original N -dimensional, generallynonlinear behavior to be studied as a 1-dimensional ordinary differential equation. In exchangefor this tremendous reduction in dimensionality, the magnitude of the perturbations is required tobe uniformly bounded in time by (cid:15) , where 0 < (cid:15) (cid:28) (cid:15) . Practically, the reduction (2) can be applied to make predictions about systembehavior in response to larger magnitude inputs with the understanding that its performance willbegin to degrade as the magnitude of the input becomes larger; precise bounds on the allowable U ( t ) are usually not known a priori and vary on an application-by-application basis.In an effort to extend the applicability of phase reduction when larger magnitude inputs areconsidered, a variety of reduction algorithms have been proposed that take into account the transientdynamics in directions transverse to the periodic orbit. For example, the notions of entrainmentmaps [11], operational phase coordinates [59], functional phase response curves [9], local orthogonalrectification [29], stochastic phase [5], and average isophase [47] rely on different definitions of‘phase’, with each being well-suited for specific applications. Other recent works have attempted touse phase-amplitude coordinate systems to identify reduced order equations that are valid beyondfirst order accuracy in (cid:15) , [58], [16], [57], [44], [61]. While such techniques generally provide betterresults than comparable first order accurate reduction methods, they still suffer from limitationsthat preclude the consideration of medium-to-large magnitude inputs. Koopman analysis is anemerging reduction framework that can be used to represent a nonlinear dynamical oscillator (ormore generally any nonlinear dynamical system) with a linear, but infinite dimensional operator[7], [33], [22]. However, it can be difficult to identify a suitable finite dimensional basis that fullycaptures the dynamics of the Koopman operator.Recently, an adaptive phase-amplitude reduction strategy was proposed in [54] that uses thestandard definition of asymptotic phase based on isochrons [64], [18] and also considers the slowlydecaying isostable coordinates [63] which represent level sets of the slowest decaying Koopmaneigenfunctions [32]. By considering a family of limit cycles that emerge for different parameter sets,the nominal parameter set can be adaptively chosen in order to limit the error in the reduced orderequations. As shown in [54], provided the Floquet exponents of the truncated isostable coordinatesare O (1 /(cid:15) ) terms, the adaptive phase-amplitude reduction is valid to O ( (cid:15) ) accuracy even when theinput U ( t ) is an O (1) term.Because the adaptive phase-amplitude reduction framework allows for the consideration of par-ticularly large magnitude inputs, this reduction framework is an attractive candidate for controlof oscillation timing in situations where the application of large magnitude inputs is unavoidable.In this work, optimal control of oscillation timing in the presence of a large magnitude entrainingstimulus will be considered. This situation can be used to represent circadian rhythms that areentrained to a 24-hour light-dark cycle. In humans, a population of roughly 10,000 coupled neu-rons referred to as the suprachiasmatic nucleus is responsible for maintaining circadian time [36],[43]. The body’s circadian cycle is nominally entrained to a 24-hour light-dark cycle in order toregulate sleep and wake timing among other physiological processes. Circadian misalignment, morecommonly referred to as jet-lag, represents a disruption to this nominal steady entrainment andresults from a mismatch between the environmental time and one’s own circadian clock [46].Numerical models of circadian oscillations have been helpful for developing post-travel treatmentprotocols to aid in recovery from circadian misalignment [3], [48], [10], however, optimal controlalgorithms can typically only be implemented on in simplified, low-dimensional models. For more2omplicated and high-dimensional models, standard phase reduction techniques typically fail. Thereare two fundamental reasons for this: the first is that circadian rhythms result from the collectiveoscillation of a large number of coupled oscillators. Phase reduction methods can be used for thesepopulation oscillations using techniques described in [25], [30], [23] and [15], however, the inter-oscillator coupling usually results in amplitude coordinates that decay slowly. Practically, theseamplitude coordinates represent shifts in the distribution of phases of the individual oscillators;without the ability to accurately account for these time-varying phase distributions only particularlysmall magnitude inputs can be considered. The second difficulty that arises when using phasereduction methods for applications involving circadian oscillations is that the nominal 24-hourlight-dark cycle itself is generally be considered a strong perturbation that shifts the oscillator statefar from its nominal limit cycle [12]. For all but a small subset of oscillator models with dynamicsthat can be greatly simplified when certain assumptions are made [27], [39], it is difficult to predicthow the aggregate behavior will respond to the large magnitude entraining stimuli precluding theconsideration of control inputs to manipulate the phase. For these reasons, reduction techniquesthat allow for large inputs are essential for applications that involve circadian oscillations.There are two primary focal points of this work. The first is the development and assessment ofstrategies to implement the adaptive phase-amplitude reduction in an optimal control frameworkto manipulate oscillation timing and speed recovery from circadian misalignment. In the results tofollow, by making appropriate assumptions about the characteristics of the adaptive reduction, acalculus of variations approach can be used to identify optimal control inputs in the presence of alarge magnitude entraining stimulus. This strategy can be used to identify inputs that engendersignificantly larger time shifts than other recently proposed oscillation timing control frameworks.The secondary focus of this work is the development of data-driven techniques for identification ofthe necessary terms of the adaptive phase-amplitude reduction that are valid in situations wherenumerical models are either unreliable or unavailable. By leveraging techniques developed in [56]a general strategy is proposed to accurately determine necessary terms of the adaptive reductionsolely from observable data.The organization of this paper is as follows: Section 2 gives necessary background informationabout the adaptive phase-amplitude reduction method used to identify reduced order models inthis work. Section 3 considers an optimal control formulation using the adaptive phase-amplitudereduction for manipulation of oscillation timing in the presence of an external entraining stimulus.Analogous, previously proposed formulations are also considered that use phase-only and firstorder accurate phase-amplitude reductions in order to provide comparisons with the proposedoptimal control formulation. Section 4 provides results for the optimal control framework whenusing a relatively simple, three-dimensional model of circadian oscillations. Section 5 considersa population-level model of coupled oscillators that give rise to a collective oscillation. Here, astrategy is presented for identifying the necessary terms of the adaptive reduction from a singlesystem observable and the optimal control strategy is subsequently implement to identify efficientstimuli to hasten recovery from circadian misalignment. Section 6 provides concluding remarks. In the analysis to follow, a phase reduction (2) will be considered that explicitly incorporates a setof nominal parameters. In this context, consider an ordinary differential equation of the form˙ x = F ( x, p ) + U ( t ) , (3)where x , F , and U are defined identically to the terms from (1) and p ∈ R M is a set of nominalparameters. Suppose that for a constant choice of p , (3) has a stable periodic orbit x γp . Provided3hat U ( t ) is sufficiently small, phase reduction is a well established technique [14], [27], [64] thatcan be used to analyze the behavior of (3) in the weakly perturbed limit according to˙ θ = ω ( p ) + Z ( θ, p ) · U ( t ) , (4)where ω and Z are defined identically to the terms from (2). In situations where larger magnitude inputs must be considered, phase reduction can still be usedbut one must also consider the transient dynamics in directions transverse to the limit cycle. Variousphase-amplitude coordinate frameworks have been developed for this purpose [58], [49], [53], [29],[5], [59], [44]; this work will focus on strategies that use isostable coordinates to characterize theamplitude coordinates [32], [58], which represent level sets of the slowest decaying eigenfunctionsof the Koopman operator [32]. Each isostable coordinate represents the magnitude of a particularKoopman eigenmode. Intuitively, states that correspond to larger magnitude isostable coordinateswill take longer to decay to the periodic orbit. For the slowest decaying Koopman eigenfunctions,explicit definitions can be given for the isostable coordinates in the basin of attraction of a fixed point[32] or periodic orbit [58], however, faster decaying isostable coordinates must be defined implicitlyaccording to their respective Koopman eigenfunctions with decay rates that are governed by theirassociated Floquet exponents. Previous work [63] (cf. [8]) investigated strategies for computing thebehavior of the isostable coordinates for models of the form (3) in response to perturbations in aneighborhood of the limit cycle yielding reduced order equations of the form˙ θ = ω ( p ) + Z ( θ, p ) · U ( t ) , ˙ ψ j = κ j ( p ) ψ j + I j ( θ, p ) · U ( t ) ,x ( θ, ψ , . . . , ψ β ) = x γp ( θ ) + β (cid:88) k =1 g k ( θ ) ψ j ,j = 1 , . . . , β, (5)where ψ j is the j th isostable coordinate with corresponding Floquet exponent κ j , g k ( θ ) are as-sociated Floquet eigenfunctions, and I j ( θ, p ) is the isostable response curve that represents thegradient of the isostable coordinates evaluated on the periodic orbit. In (5), the rapidly decayingamplitude coordinates are generally truncated so that only β ≤ N − Z , and each I j and g k can be performed using the ‘adjoint method’ [6] and relatedequations described in [57]. The so-called ‘direct method’ [21], [38] and related data-driven tech-niques [60], [56] have been developed for inferring the necessary terms of (5) from experimentaldata when the underlying model equations are unknown.Much like (4), Equation (5) is only valid provided the state remains sufficiently close to theperiodic orbit. Unlike (4) alone, however, (5) can be used to provide information about the tran-sient behaviors that characterize the amplitude dynamics. It is generally assumed that under theapplication of input with O ( (cid:15) ) magnitude that each ψ j remains an O ( (cid:15) ). In this case, (5) is valid toleading order (cid:15) . Recent work has investigated isostable reduced equations that are valid to secondorder accuracy [58], [60] and beyond [57], however, these reduction frameworks still require themagnitude of the applied inputs to be sufficiently small. Other reduction frameworks have been de-veloped that are valid for inputs with arbitrary magnitude provided that they are sufficiently slowly4arying [28], [40] or rapidly varying [42], [52]. Nevertheless, few general reduction frameworks areavailable that are valid for arbitrary, large magnitude inputs. Isostable reduced equations of the form (5) typically assume that the nominal system parametersare constant in the analysis of a perturbed limit cycle. Recent work [54] has leveraged the phase-isostable coordinate reduction framework (5) to develop an adaptive reduction that actively setsthe nominal system parameters in an effort to keep isostable coordinates low, thereby allowingfor inputs without O ( (cid:15) ) constraints (either in magnitude or rate of change). As discussed in [54],to implement an adaptive phase-isostable reduction, first suppose that in some allowable range ofnominal system parameters p ∈ R M that x γp ( t ) is continuously differentiable with respect to both t and p . One can then rewrite Equation (3) as˙ x = F ( x, p ) + U e ( t, p, x ) , (6)where U e ( t, p, x ) = U ( t ) + F ( x, p ) − F ( x, p ) . (7)Intuitively, F ( x, p ) represents the underlying dynamics for a given choice of system parameters p , and U e represents the effective input. Allowing p to be variable, and assuming that θ ( x, p ) and ψ j ( x, p ) are continuously differentiable in a neighborhood of the limit cycle for all x and p , one cantransform (6) to phase and isostable coordinates via the chain rule dθdt = ∂θ∂x · dxdt + ∂θ∂p · dpdt ,dψ j dt = ∂ψ j ∂x · dxdt + ∂ψ j ∂p · dpdt , (8)for j = 1 , . . . , β . Above, provided that the state is close to the periodic orbit x γp , ∂θ∂x · dxdt = ω ( p ) + Z ( θ, p ) · U ( t ) and ∂ψ j ∂x · dxdt = κ j ( p ) ψ j + I j ( θ, p ) · U ( t ) as given by Equation (5). For theremaining terms, as explained in [54], one can write D ( θ, p ) ≡ ∂θ∂p = (cid:104) Z ( θ, p ) · ∂ ∆ x p ∂p . . . Z ( θ, p ) · ∂ ∆ x p ∂p M (cid:105) T ,Q j ( θ, p ) ≡ ∂ψ j ∂p = (cid:104) I j ( θ, p ) · ∂ ∆ x p ∂p . . . I j ( θ, p ) · ∂ ∆ x p ∂p M (cid:105) T , (9)where ∂ ∆ x p ∂p i = lim dp i → x γp − x γp + dp i dp i . (10)Note that Equation (10) represents the change in the nominal periodic orbit in response to a changein the parameter p i . In Equation (9) and (10), all derivatives are evaluated at p and θ on the limitcycle x γp . Substituting both (9) and the phase and isostable dynamics from (5) into (8) yields theadaptive phase-isostable reduction˙ θ = ω ( p ) + Z ( θ, p ) · U e ( t, p, x ) + D ( θ, p ) · ˙ p, ˙ ψ j = κ j ( p ) ψ j + I j ( θ, p ) · U e ( t, p, x ) + Q j ( θ, p ) · ˙ p,j = 1 , . . . , β, ˙ p = G p ( p, θ, ψ , . . . , ψ β ) , (11)5here the function G p is designed to actively choose p (and by extension ˙ p ) in a manner that keepsthe isostable coordinates small. As explained in [54], provided G p can be designed so that ψ j remain O ( (cid:15) ) for j ≤ β and that the neglected isostable coordinates have sufficiently large magnitude Floquetexponents, Equation (11) is accurate to O ( (cid:15) ) provided that U e is an O (1) term. Recalling that thestandard phase-isostable reduction (5) assumes the input is an O ( (cid:15) ) term, the adaptive isostablereduction represents a substantial improvement that allows for significantly larger magnitude inputsto be considered. It is worth emphasizing that the underlying model (3) does not need to havetime-varying parameters in order for the adaptive reduction (11) to be implemented. Rather, theadaptive reduction actively changes the nominal parameter set p so that the state is close to theperiodic orbit x γp thereby limiting the magnitude of the isostable coordinates. Both phase reduction (4) and phase-amplitude reduction have been applied fruitfully to applicationsthat involve control and analysis of behaviors that emerge in limit cycle oscillators [50], [1], [61],[66], [41]. The effective reduction in dimension that results from phase reduction can allow for theapplication of control and analysis techniques that would otherwise be intractable. However, theresults are only applicable when sufficiently weak perturbations are applied limiting the practicalutility in many applications. Here, a prototype problem is considered for identifying an optimalcontrol input to advance or delay the oscillation timing in an oscillator subject to an externalentraining stimulus. This problem is motivated by the search for jet-lag mitigation protocols thatcan be used to limit the negative effects of circadian misalignment by allowing one’s circadian cycleto reentrain rapidly to a new time zone [3], [48], [10]. Related problems were considered in [34]and [37] using standard phase reduction methods, and in [35] (resp. [58]) using first (resp. second)order accurate phase-amplitude reduced equations. As shown in the results to follow, the problemformulation using the adaptive phase-amplitude reduced equations allows for inputs that are farlarger in magnitude than those considered previously without sacraficing accuracy – consequently,inputs that provide significantly larger magnitude changes in oscillation timing can be obtained.To begin, consider the adaptive phase-isostable reduced equations from (11) with only a singleisostable coordinate ˙ θ = ω ( p ) + Z ( θ, p ) · U e ( t, p, x ) + D ( θ, p ) ˙ p, ˙ ψ = κ ( p ) ψ + I ( θ, p ) · U e ( t, p, x ) + Q ( θ, p ) ˙ p, ˙ p = G p ( p, θ, ψ ) . (12)Above, for convenience of notation, the subscripts on the terms involving isostable coordinates havebeen omitted because there is only one isostable coordinate. Additionally, it will be assumed that p ∈ R . In the analysis to follow, it will be assumed that the input can be written as U ( t ) = δu ( t ),where δ ∈ R N is a constant vector and u ( t ) ∈ R is a time-varying control signal. In other words, U ( t ) a rank-1 perturbation. It will also be assumed that F ( x, p ) = F ( x ) + pδ . This situationresults, for instance, when allowing the input to also be considered as a time-varying parameter inthe adaptive reduction. Under these assumptions, Equation (12) simplifies to˙ θ = ω ( p ) + z ( θ, p )( u ( t ) + p − p ) + D ( θ, p ) ˙ p, ˙ ψ = κ ( p ) ψ + i ( θ, p )( u ( t ) + p − p ) + Q ( θ, p ) ˙ p, ˙ p = G p ( p, θ, ψ ) , (13)6here z ( θ, p ) = Z ( θ, p ) · δ and i ( θ, p ) = I ( θ, p ) · δ . To further simplify (13), it will be explic-itly assumed that | Q ( θ, p ) | > ν for all allowable θ and p with ν >
0. In this situation, taking G p ( p, θ, ψ ) = − R ( θ, p )( u + p − p ) where R ( θ, p ) = i ( θ,p ) Q ( θ,p ) , the isostable dynamics from equation(13) simplify to ˙ ψ = κ ( p ) ψ which has a single globally stable equilibrium at ψ = 0. Noting that G p does not depend on ψ , the isostable coordinate dynamics can be ignored from (13) yielding˙ θ = ω ( p ) + z ( θ, p )( u ( t ) − p ) − D ( θ, p ) R ( θ, p )( u − p ) , ˙ p = − R ( θ, p )( u − p ) , (14)where p = 0 is assumed for simplicity of exposition. To formulate the optimal control problem,consider a T ext -periodic external input u ( t ) = u nom ( t s ) applied to Equation (3) where t s = (cid:40) t, for t ≤ t ,t + ∆ t, for t > t . (15)Suppose that for t ≤ t , the oscillator is fully entrained to the periodic orbit so that x ( t ) = x ( t + T ext )for all t + T ext <
0. This fully entrained orbit will be denoted by x γ ent ( t s ). In the context ofcircadian oscillations, t s would represent the environmental time that controls a 24-hour light-darkcycle u nom ( t s ). Suppose that at t = t , the environmental time instantaneously shifts by ∆ t timeunits, for instance, representing a sudden shift across multiple time zones. The control objectiveis to identify a control input u ( t ) = u nom ( t + ∆ t ) + ∆ u ( t ) so that x ( t + T f ) = x γp , ent ( t + ∆ t )that minimizes the cost functional C = (cid:82) t + T f t ∆ u ( t ) dt subject to the constraints on the allowablecontrol ∆ u min ≤ u ( t ) + u nom ( t + ∆ t ) ≤ ∆ u max . Intuitively, this control problem seeks to identifyan external input ∆ u ( t ) so that the system is fully entrained after T f time units to the time-shiftedentraining stimulus.Working in the adaptive reduction framework from (14) this control problem can be approachedby first defining the Hamiltonian associated with the cost functional H (Φ , ∆ u, Λ , t ) = ∆ u ( t ) + λ (cid:2) − R ( θ, p )(∆ u ( t ) + u nom ( t + ∆ t ) − p ) (cid:3) + λ (cid:2) ω ( p ) + z ( θ, p )(∆ u ( t ) + u nom ( t + ∆ t ) − p ) − D ( θ, p ) R ( θ, p )(∆ u ( t ) + u nom ( t + ∆ t ) − p ) (cid:3) , (16)where Φ ≡ (cid:2) θ p (cid:3) T contains the state variables and Λ ≡ (cid:2) λ λ (cid:3) T represents Lagrange multipliersthat force the dynamics to satisfy the dynamics of the adaptive reduction. According to Pon-tryagin’s minimum principle [24], the control that minimizes the cost functional C will minimizethe Hamiltonian for all admissible ∆ u ( t ) with the dynamics of the state variables and Lagrangemultipliers subject to ˙ x = ∂ H ∂ Λ , (17)˙Λ = − ∂ H ∂x . (18)Evaluation of (17) returns the state equations of the adaptive reduction from (14). Evaluation of(18) yields˙ λ = − λ (cid:2) z θ (∆ u ( t ) + u nom ( t + ∆ t ) − p ) − ( D θ R + DR θ )(∆ u ( t ) + u nom ( t + ∆ t ) − p ) (cid:3) + λ R θ (∆ u ( t ) + u nom ( t + ∆ t ) − p ) , ˙ λ = − λ (cid:2) ω p + z p (∆ u ( t ) + u nom ( t + ∆ t ) − p ) − z − ( D p R + DR p )(∆ u ( t ) + u nom ( t + ∆ t ) − p ) + DR (cid:3) − λ (cid:2) R − R p (∆ u ( t ) + u nom ( t + ∆ t ) − p ) (cid:3) , (19)7here, for instance, the notation z θ corresponds to the partial of z with respect to θ evaluated atboth θ and p . Noting that the Hamiltonian (16) is quadratic in ∆ u , the admissible control ∆ u ( t )that minimizes the Hamiltonian at any given moment is∆ u = min (cid:18) max (cid:18) − z ( θ, p ) λ + λ D ( θ, p ) R ( θ, p ) + R ( θ, p ) λ , ∆ u min (cid:19) , ∆ u max (cid:19) . (20)In the adaptive reduction framework, letting θ ent ( t ) and p ent ( t ) correspond to the θ and p valueson the entrained periodic orbit, the boundary conditions prescribed by the control problem are θ ( t ) = θ ent ( t ), θ ( t + T f ) = θ ent ( t + ∆ t ), p ( t ) = p ent ( t ), and p ( t + T f ) = p ent ( t + ∆ t ). Inorder to solve this two-point boundary value problem, it is necessary to identify the correct choiceof λ ( t ) and λ ( t ) that yield the prescribed final conditions. This is accomplished in this work byproviding an initial guess and subsequently updating the initial values of these Lagrange multipliersusing a Newton iteration until convergence is achieved. In addition to the optimal control formulation that uses the adaptive phase-amplitude reduction,two additional control formulations used for comparison purposes. The first control strategy usesthe first order accurate phase-isostable reduced equations from (5) with dynamics that follow˙ θ = ω + z ( θ )(∆ u ( t ) + u nom ( t s )) , ˙ ψ = κψ + i ( θ )(∆ u ( t ) + u nom ( t s )) , (21)where ω and κ are evaluated at p , U ( t ) = δu ( t ) as defined directly above (13), with z ( θ ) = Z ( θ, p ) · δ and i ( θ ) = I ( θ, p ) · δ denoting the effective phase and isostable response curves. Compared withthe adaptive reduction (14), Equation (21) does not allow for adjustments of the underlying systemparameter p , but rather, explicitly takes into account the isostable coordinates defined in referenceto the periodic orbit x γp . Such a formulation allows for the definition of a cost functional, similarto the one proposed in [35], that balances the tradeoff between the magnitude of the isostablecoordinate and the control input C = (cid:90) t + T f t ∆ u ( t ) + βψ ( t ) dt, (22)where β is a positive constant. Using the same control objective, same limits on the optimal control,and same initial and final conditions as those from the previous sections, one can follow an identicalset of arguments to write the Hamiltonian associated with the cost functional (22) H (Φ , ∆ u, Λ , t ) = ∆ u ( t ) + βψ + λ (cid:2) ω + z ( θ )(∆ u ( t ) + u nom ( t + ∆ t )) (cid:3) + λ (cid:2) κψ + i ( θ )(∆ u ( t ) + u nom ( t + ∆ t )) (cid:3) , (23)where Φ = (cid:2) θ ψ (cid:3) T . Associated optimal trajectories must satisfy (21) along with˙ λ = − λ z θ (∆ u ( t ) + u nom ( t + ∆ t )) − λ i θ (∆ u ( t ) + u nom ( t + ∆ t )) , ˙ λ = − λ κ − βψ, ∆ u = min (cid:18) max (cid:18) − λ z ( θ ) − λ i ( θ )2 , ∆ u min (cid:19) , ∆ u max (cid:19) . (24)8n the phase-isostable reduced framework, letting θ ent ( t ) and ψ ent ( t ) correspond to the θ and p values on the entrained periodic orbit, the boundary conditions prescribed by the control problemare θ ( t ) = θ ent ( t ), θ ( t + T f ) = θ ent ( t + ∆ t ), ψ ( t ) = ψ ent ( t ), and ψ ( t + T f ) = ψ ent ( t + ∆ t ).The second alternative control formulation considered only uses the phase dynamics, completelyignoring all amplitude coordinates. Such a control strategy only considers the top equation from(21) with a cost function C = (cid:90) t + T f t ∆ u ( t ) dt. (25)The corresponding Hamiltonian is H ( θ, ∆ u, λ , t ) = ∆ u ( t ) + λ (cid:2) ω + z ( θ )(∆ u ( t ) + u nom ( t + ∆ t )) (cid:3) , (26)with associated optimal control trajectories that must satisfy˙ θ = ω + z ( θ )(∆ u ( t ) + u nom ( t + ∆ t )) , ˙ λ = − λ z θ ( θ )(∆ u ( t ) + u nom ( t + ∆ t )) , ∆ u = min (max ( − λ z ( θ ) / , ∆ u min ) , ∆ u max ) . (27)Without information about the isostable dynamics, the boundary conditions only involve the phasecoordinates and are taken to be θ ( t ) = θ ent ( t ), θ ( t + T f ) = θ ent ( t + ∆ t ). Related optimal controlproblems that only account for the phase dynamics were considered in [34] and [37]. Results are presented below that use adaptive phase-isostable reduction framework in conjunctionwith the prototype optimal control problem. In the context of circadian cycles, the optimal controlproblem can be viewed as a strategy that seeks to mitigate circadian misalignment (i.e., jet-lag)caused by a sudden shift in the environmental time occurring as a result of rapid travel throughmultiple time zones. In this section, a simple three-dimensional model for circadian oscillations isconsidered. In the section to follow, a collective oscillation that represents the aggregate behaviorof 3000 coupled oscillators will be considered.For this example, consider a simple three-dimensional model for circadian oscillations [17]˙ B = v K n K n + D n − v BK + B + L c + L nom ( t s ) + ∆ L ( t ) , ˙ C = k B − v CK + C , ˙ D = k C − v DK + D . (28)In the above equations, B , C , and D , represent concentrations of mRNA of the clock gene, associ-ated protein, and nuclear form of the protein, respectively. The term L nom ( t s ) represents a nominal24-hour light-dark cycle taken to be L nom ( t s ) = L (cid:20)
11 + exp( − t s − −
11 + exp( − t s − (cid:21) , (29)9here L is the maximum uncontrolled light intensity and t s = mod( t + ∆ t,
24) with ∆ t being atime shift. The term ∆ L ( t ) is taken to be a control input with the overall light intensity boundedby L min ≤ L nom ( t s ) + ∆ L ( t ) ≤ L max . Matching the formulation given in (13), U ( t ) = δ (∆ L ( t ) + u nom ( t )) where δ = (cid:2) (cid:3) T . Remaining parameters are taken to be n = 6, v = 0 . , v =0 . , v = 0 . , v = 0 . , K = 1 , K = 1 , K = 1 , K = 1 , k = 0 .
7, and k = 0 .
7. With thisparameter set, when L c = L nom = ∆ L = 0, (i.e., in the absence of light) the resulting stable limitcycle has a period of 24.2 hours.Letting L c be the time-varying parameter used in the adaptive reduction (13), the traces of B ( θ ) on the resulting periodic orbits x γL c are shown in panel A of Figure 1. The periodic orbits areappropriately shifted so that θ = 0 corresponds to the moment that B ( θ ) reaches its peak value.For all values of L c considered, the periodic orbit has only one non-negligible Floquet exponent (theother Floquet exponent is negative and large in magnitude). Along with this Floquet exponent,the nominal period is shown as a function of L c in panel B . The resulting phase and isostableresponse curves for various values of L c are computed using an adjoint method described in [55]and shown in panels C and D, respectively. Finally, D and Q are computed according to (9) andshown in panels E and F, respectively.Figure 1: Traces of B γ ( θ ) are shown for (28) using various constant values of L c . The resultingnatural frequencies and Floquet exponents are shown in panel B. The dashed line at κ = 0 is shownfor reference, emphasizing that all considered orbits are stable. Phase and isostable response curvesused in the adaptive reduction are shown in panels C and D, respectively. The terms D and Q characterize how changing the nominal parameter influences the phase and isostable coordinateand are shown in panels E and F, respectively. Note that both the magnitude of oscillations as wellas the magnitude of the Floquet exponents become smaller as L c increases resulting less robustoscillations.For the moment, the optimal control problem will be considered taking L = 0. Additionally,it will be assumed that L max = + ∞ , and L min = −∞ . This situation represents a situation whereno entraining stimulus is applied and no bounds on the magnitude of light are present. While thissituation is not realistic, it provides insight into the limitations of each optimal control formulation.10ote that when L = 0, the system is technically not entrained to any external stimulus. Never-theless, the optimal control frameworks proposed in Section 3 can still be implemented by defining x γ ent ≡ x γp with x γp | θ =0 = x γ ent | t =0 .For various choices of ∆ t , an optimal control is computed using the adaptive reduction frame-work (with Equations (14), (19) and (20)), the first order accurate phase-amplitude reductionstrategy (with Equations (21) and (24)), and the phase-only reduction strategy (with Equation(27)). For the phase-amplitude and phase-only reduction strategies, L c is taken to be 0. Theconstant β (which sets the penalty for large isostable coordinates) is taken to be 10 − for thephase-amplitude control strategy. T f is taken to be 24 hours for all control strategies. Optimalinputs ∆ L ( t ) are computed using each reduction strategy, with the resulting inputs applied to thefull model (28).Figure 2: Optimal control results when no external entraining stimulus is applied and no con-straints are placed on ∆ L ( t ). The resulting optimal controls computed using the adaptive reduction,phase-isostable reduction, and the phase-only reduction are applied to the full model equations (28).The specified ∆ t is plotted against the resulting time change, ∆ t act in Panel A. The error, definedto be ∆ t act − ∆ t is shown in panel B. Specific inputs are shown for ∆ t = +9 and -9 hours in panelsC and D using the specified method. Resulting traces of B ( t ) when these inputs are applied tothe full model are shown in Panels E and F, respectively. The black dashed line shows the fullyentrained solution after the time shift is applied. While all methods perform reasonably well forfor | ∆ t | < t .Resulting time shifts are shown in panel A of Figure 2. Note that because there is no externalentraining stimulus in this example, mandating a time shift ∆ t is equivalent to mandating a phaseshift ∆ θ = 2 π ∆ t/ .
2, and t act is the actual resulting phase shift when the input is applied tothe full model (28). ∆ t act is computed by considering the infinite-time behavior after the input isapplied. In panel B, Error = ∆ t act − ∆ t is shown. For small values of ∆ t , all reduction frameworksproduce stimuli that achieve the control objective. As the required ∆ t shift increases, errors tend toincrease for the phase-only and the phase-amplitude control strategies. By contrast, the adaptivereduction strategy yields results that are nearly perfect, with the maximum value of | ∆ t act − ∆ t | being less than 0.4 hours for any choice of ∆ t . Panel C (resp., D) show resulting control inputs11or each strategy that result when ∆ t = +9 hours (resp. -9 hours). Panels panel E and F showcorresponding traces of B ( t ) that result when those inputs are applied to the full model (28).Next, a situation where L = 0 .
01 is considered next so that the limit cycle is nominallyentrained to a the external 24-hour light-dark cycle (29). Noting that the nominal phase responsecurves from panel C of Figure 1 can take values that are on the order of approximately 20, thisinput is quite large for this system. Such situations have traditionally been difficult to approachwith standard phase reduction techniques (as discussed in [11]) because the entraining input isoften strong enough to drive the state far from its unperturbed limit cycle, thereby invalidating theassumptions of the phase reduction. Alternative techniques that consider the dynamics of Poincar´emaps [11], [12] or those that consider the entrained orbit itself to be the nominal periodic orbit[56] have been useful in some situations. As shown in the results to follow, the control strategythat uses adaptive phase-amplitude reduction yields accurate results while the phase-only and firstorder accurate phase-amplitude reduction strategies fail.For this example, L min = 0 is chosen in optimal control problem along with L max = + ∞ tomandate L nom +∆ L ≥ t , an optimal control is computed using the adaptive reduction framework(from Equations (14), (19) and (20)), the first order accurate phase-amplitude reduction strategy(from Equations (21) and (24)), and the phase-only reduction strategy (from Equation (27)). Forthe phase-only reduction strategy β is set to zero, since the nominal 24-hour light-dark cycle makesit difficult to limit the magnitude of the isostable coordinate. In all optimal control computations, t = 0 which corresponds to the middle of subjective night. T f is taken to be 24 hours for thephase-only and adaptive reduction strategies. T f is taken to be 72 hours for for the phase-amplitudereduction strategy, it is not possible to find numerical solutions to the optimal control problem forsmaller values of T f . Resulting inputs are then applied to the full model equations (28) with resultsshown in Figure 3. Because the model is entrained to L nom ( t ) in the absence of additional input,the resulting recovery time, defined to be the time it takes for the phase to return to and staywithin one hour if its steady state behavior, is shown in Panel A. Results are also compared to thenominal, uncontrolled recovery time (i.e. taking ∆ L = 0). Note that the phase is only measuredat θ = 0, which can be observed once per cycle when B ( t ) reaches a local maximum. The recoverytime is inferred by computing the difference in time between the moment θ reaches 0 for the fullyentrained reference and the reentraining simulation and interpolating the time difference for allvalues in between. If the time difference is less than one hour the first moment θ = 0 is reachedafter input is applied for the reentraining system, the recovery time is taken to be equal to T f .The optimal controls computed according to the phase-isostable and the phase-only reductionstrategies generally have only a small influence on the recovery times and even occasionally result inworse outcomes than if no input was applied at all. Conversely, the when the control identified whenusing the adaptive phase-amplitude reduction is applied to the full model equations, the recoverytime exactly 24 hours indicating the system is recovered by the time the first measurement of θ = 0is made after the optimal stimulus is applied. Panel B of Figure 3 shows L c when the optimalcontrol is applied to the adaptive reduction. Panels C and D (resp., E and F) show the optimalcontrol inputs and traces of B ( t ) for a time shift of ∆ t = +12 hours (resp. -8 hours). Black linesalso show traces of B ( t ) and the nominal light-dark cycle input during the uncontrolled recovery.12igure 3: Optimal Control results with a nominal external entraining stimulus and constraints onthe applied stimulus are considered. Panel A shows recovery times when applying the resultingoptimal inputs determined from each of the three optimal control formulations to the full modelequations (28). The uncontrolled recovery times are also showed for reference. The formulationwith the adaptive reduction significantly outperforms the other two control formulations. Indeed,the phase-only and first order accurate phase-isostable formulations often yield control inputs that increase the recovery time. Panel B shows the value of p , i.e., the parameter that corresponds tothe nominal periodic orbit during the application of the optimal stimulus for various time shifts ∆ t .Notice that p increases to much larger values for all but the smallest magnitude choices of ∆ t . Thesechanges in p represent large deviations from the nominal limit cycle which has p = 0. The adaptivereduction has the ability to explicitly incorporate these shifts, however, the phase-amplitude andphase-isostable reductions cannot explicitly account for these shifts. This fundamental differencebetween the reduction frameworks leads to the difference in efficacy between the resulting optimalinputs. Panels C and D show the recovery for the optimal controls identified by each of the optimalcontrol formulations for ∆ t = +12 hours. The black curve also shows the nominal recovery when∆ L ( t ) = 0. The dashed line in panel C represents the limit cycle that is fully entrained to theshifted nominal input L nom ( t + ∆ t ). Panels E and F show the same information for a time shift of-8 hours. Here, a model for the behavior of a population coupled circadian oscillators will be considered:˙ B i = v K n K n + D ni − v B i K + B i + h c KFK c + KF + α i ( L c + L nom ( t s ) + ∆ L ( t )) , ˙ C i = k B i − v C i K + C i , ˙ D i = k C i − v D i K + D i , ˙ E i = k B i − v E i K + E i , i = 1 , . . . , N. (30)13he above model considers a population of N = 3000 coupled oscillators of the same form as (28)with the addition of the dynamics of a neurotransmitter E i that allows communication betweencells. Assuming that the spatial transmission of the neurotransmitter is fast relative to the 24-hour scale of oscillation, effective mean-field coupling is assumed with F ≡ (1 /N ) (cid:80) Nk =1 E i , whichenters into the equation for the variable B . Additionally, in (30) each neuron has an intrinsicsensitivity to light, α i , drawn from the distribution α i = max(1 + 0 . N (0 , ,
0) where N (0 , n = 5 , v . , v . , k = 0 . , v = 0 . , k = 0 . , v = 0 . , k =0 . , v = 1 , K = 1 , K = 1 , K = 1 , K = 1 , K = 1 , h c = 1 , K c = 1, and K = 0 .
5. In order toincorporate heterogeneity in the model, the parameters k , k , k , v , v , and v are drawn from anormal distribution with the mean being the nominal parameter value and a standard deviationequal to 0.01. In all simulations, independent and identically distributed zero mean white noisewith intensity 0.0002 is added to the variable B i for each oscillator. Like in the single oscillatorequations (28), the external light-dark cycle L nom ( t s ) in (30) takes the form (29).While each of the necessary terms of the reduction in (13) could be computed numerically, astrategy for inferring the required terms from output data will be considered here, as would be thecase for an experimental system. The control strategy from Section 3 will then be applied to theresulting adaptive reduction. When the underlying dynamical equations are known, it is relatively straightforward to compute thenecessary terms of the adaptive reduction (13) using the ‘adjoint method’ [6] and related equationsfrom [57] as described in Section 2. However, when the the dynamical equations are unknown,the terms of the adaptive reduction must inferred from data. The ‘direct method’ [21], [38] is awell-established technique for identifying the phase response curve from (4). Previous work [56]developed a data-driven technique for computing the terms of the phase-amplitude reduction from(5). This strategy will be expanded here to compute the necessary terms of the adaptive reductionfrom (13). Here, it is assumed that the input also the sole time-varying parameter. While (13) hasonly one isostable coordinate, this technique is relatively straightforward to extend to situationswhere multiple isostable coordinates are involved.For concreteness, the model (30) will be used assuming light perturbations ∆ L ( t ) can be givenand that the single output F ( t ) (i.e., the mean value of E i for the population) can be measured.For the given parameter set it will be assumed that Equation (30) has only one dominant isostablecoordinate. Additionally, L c will be taken as the time-varying paramater. In this situation, it ispossible to write the equations for the adaptive phase-amplitude reduced model in the form (13).The terms of the adaptive reduction will be inferred using the procedure detailed below. Step 1:
With a static choice of L c , and taking L nom ( t s ) + ∆ L ( t ) = 0 for all time, the procedureintroduced in [56] can be used to identify ω ( L c ), κ ( L c ). This strategy can also be used to identify ψ ( t ) and θ ( t ) provided θ ( t ) ≈ F ( t ) to F ( t + T ).This strategy is summarized in Figure 4. To apply this strategy, the model is first simulated untiltransient behaviors have died out. A threshold is chosen to denote θ ≈ θ ≈ F ( t ) crosses 0.045 with a positive slope). Over multiple oscillations,the average period at which this threshold is crossed is taken to be the period of oscillation, T . Aset of delay embeddings, each which start the moment that F ( t ) crosses 0.045 and end T hourslater, is recorded and the average value of the output F over these measurements is used as anapproximation for the stable periodic orbit F γ ( t ) as shown in panel B of Figure 4. Once the periodic14rbit is obtained, the recovery to the periodic orbit from a perturbed initial condition is consideredin order to determine the isostable coordinates. As shown in panel A, the coupling strength K isdecreased 50 percent for t ∈ [0 , t > K is set back to its nominal value, and aseries of delayed embeddings from t j to t j + T are taken, where each t j is chosen so correspond toa moment that F ( t j ) crosses 0.045 with a positive slope. This procedure is repeated over multipletrials. The periodic orbit is subtracted from the resulting embeddings with individual traces shownin panel C. To proceed, let y j ∈ R q correspond to the j th delay embedding from panel C. In otherwords, y j is comprised of the output from t j to t j + T after F γ is subtracted. Here, q = ( T /δt ) + 1where δt is the sampling rate.Figure 4: A data-driven strategy for identifying phase and isostable coordinates from the largecircadian model (30). In panel A, the model is perturbed from its limit cycle by decreasing K bya factor of 2 for 200 hours (red curve). The subsequent recovery takes approximately 400 hours(blue curve) and the black curve represents the behavior once the model has recovered to thestable limit cycle. As shown in panel B, the periodic orbit F γ ( t ) (black line) is taken to be theaverage of multiple cycles once transient behavior has died out. Using data obtained from the blueportion of panel A, multiple delayed embeddings are obtained and shown in panel C that representthe deviation from the periodic orbit on a particular cycle. POD is performed on the data frompanel C and identifies two characteristic modes. The POD modes are transformed to represent thedeviation from the periodic orbit in a basis of phase and isostable coordinates with modes shownin panel E. As described in the text, this procedure allows one to measure the isostable coordinatesat specific times, and a the slope of the linear regression shown in Panel F represents the Floquetmultiplier corresponding to the isostable coordinate.Proper orthogonal decomposition (POD) [4], [19], [51] is applied to the data shown in panel Cin order to identify a small number of representative modes. In this case, 2 modes are sufficientto represent the data from panel C, denoted by φ ∈ R q and φ ∈ R q and shown in panel D asblue and red curves, respectively. The transformation (cid:2) µ µ (cid:3) = (cid:2) φ φ (cid:3) A where A is a 2 × µ = (cid:2) φ φ (cid:3) † g γ where g γ ∈ R q corresponds to dF γ /dt taken at the sampling rate δt and † represents the Moore-Penrose pseudoinverse. Plots of µ and µ are shown in panel E in blue and red, and are proportional to the shifts in the outputresulting from a shift in θ and ψ , respectively. Furthermore as described in [56], the vectors η and η defined according to, A − (cid:2) φ φ (cid:3) T = (cid:2) η η (cid:3) T can be used to compute the phase and15sostable coordinates at specific instances in time according to the relationships cη T y j = θ ( t j ) ,η T y j = ψ ( t j ) , (31)where c is a constant that can be determined using a strategy discussed in [56]. This informationcan be used to track ψ over successive periods during the recovery as plotted in panel F. The slopeof a linear regression of this data gives an approximation of the Floquet multiplier, Λ, correspondingto ψ and the Floquet exponent can be computed according to the relationship κ = log(Λ) /T .Note that the procedure from Step 1 represents an implementation of the procedure given inSection IV,A from [56] for a constant choice of L c . The interested reader is referred to [56] for amore detailed description of this method. This strategy is repeated to determine the associatedterms of the phase and isostable coordinate reduction for various values of L c . Step 2:
For a given choice of L c the terms of z and i can be inferred directly by applying a seriesof pulse inputs and measuring the resulting change in the phase and isostable coordinates. Thiscan be accomplished with a strategy akin to the direct method [21], [38] whereby a pulse input ofmagnitude dL is applied for a short duration of time as shown in panel B of Figure 5. Completecycles of the output starting at θ ≈ F ( t ) = 0 .
045 threshold)beginning at t and t preceding and following the input can be extracted as shown in panel A, andthe relationships from (31) can be used to determine θ ( t ), θ ( t ), ψ ( t ), and ψ ( t ). Assuming thepulse input was applied at t = 0, the shift in phase can be computed by recalling that the phasewould have simply increased at the rate ω had the input not been applied. The deviation, ∆ θ ,from this expectation can be assumed to be the change in phase caused by the input. Likewise, byassuming that ψ decays exponentially towards zero at a rate governed by κ in the absence of input, ψ can be computed immediately before and after the applied input to infer the change in isostablecoordinate ∆ ψ . For the chosen parameter L c , z ( θ, L c ) ≈ ∆ θdLdt and i ( θ, L c ) ≈ ∆ ψdLdt , where dt is theduration of the perturbation applied at phase θ . This process can be repeated for multiple trialsto obtain approximations of i ( θ, L c ) and z ( θ, L c ) as shown in Panels C and D of Figure 5 and theresulting data can be fit using an appropriate basis. Step 3:
Both D and Q can be inferred using a strategy similar to the one described in Step 2. Inthis case, instead of a brief pulse, the total applied light is shifted from L c to L c + dL at a knowntime as shown in panel F of Figure 5. This shift can be viewed as a sudden change in the nominalparameter by dL centered at L c + dL/
2; in other words a pulse to ˙ L c of magnitude dL . Consideringthe structure of the reduced Equations (13), this pulse in ˙ L c can be used to infer both D and Q ina manner similar to the direct method for identifying z and i . Similar to the procedure used in step2, cycles beginning at t and t preceding and following the step function change in input as shownin panel E of Figure 5 can be isolated from the output. Subsequently, Equations (31) can be usedto identify θ ( t ), θ ( t ), ψ ( t ), and ψ ( t ). Note that for each trial θ ( t ) and ψ ( t ) (resp. θ ( t ) and ψ ( t )) must be found using η and η that are obtained from Step 1 taking the nominal parameterto be L c (resp. L c + dL ). Using the same reasoning employed as part of Step 2, the change inphase, ∆ θ , and isostable coordinate, ∆ ψ , caused by the parameter change can be inferred yieldingindividual data points for D ( θ, L c + dL/
2) = ∆ θ/ ∆ p and Q ( θ, L c + dL/
2) = ∆ ψ/ ∆ p . This processcan be repeated to obtain multiple data points; curves can be fit to the resulting data as shown inpanels G and H of Figure 5. Step 4:
Steps 3 and 4 can be repeated for various choices of L c to obtain a series of curves usedin the reduced equations of the form (13). Resulting information obtained when applying theprocedure from Steps 1 through 4 to the model from (30) are shown Figure 6. Numerical trials16igure 5: Inferring the phase and isostable response to pulse and step function inputs to infer theterms of the reduction (13). For a nominal value of L c , a short pulse can be applied as in panelB. Data from a complete cycle occurring both before and after the perturbation (red curves inpanel A) are used to identify the phase and isostable coordinates at t and t . This information isthen used to obtain a discrete measurement of z ( θ, L c ) and i ( θ, L c ), and this process is repeatedmultiple times. Each trial yielding a single measurement of both z and i which are representedby black dots in panels C and D. The blue lines are obtained by fitting the sinusoidal basis of theform (cid:80) n =0 [ a n sin( nθ ) + b n cos( nθ )]; the resulting fits are taken to be the response curves. Likewise,panels E through H illustrate a similar procedure for obtaining D and Q . Instead of applying pulses,step function inputs are applied (panel F), complete cycles before and after the step function areextracted (panel E) and used to infer the change in phase and isostable coordinates resulting fromthe input. Each trial gives a single data point for D and Q (black dots in panels G and H). Thisprocess is repeated for multiple trials and the same sinusoidal basis is used to fit the data (bluelines).are performed taking L c ∈ {− . , . , . , . , . } , and all other curves are obtained throughinterpolation. 17igure 6: Resulting terms of the adaptive reduction (13) resulting from the data-driven inferencestrategy from Section (5.1). For all limit cycles, θ = 0 corresponds to the moment that F ( t ) crosses0.045 with a positive slope. As L c increases, the magnitude of the nominal oscillations decrease asshown in panel A. Panels B and C give the natural frequency and Floquet exponent as a functionof L c . Panels D-F show specific terms of the adaptive reduction for various values of L c . The terms of the adaptive reduction inferred using the strategy from Section 5.1 and shown inFigure 6 are used to implement the optimal control strategy for manipulating oscillation timing asdescribed in Section 3. The nominal value of L c is taken to be zero and the nominal magnitude ofthe light-dark cycle is taken to be L = 0 .
01. Bounds on the total allowable light input are chosento be L min = 0 and L max = ∞ so that negative light inputs are not possible. Recalling that thechoice of ∆ t in the control formulation corresponds to a sudden change in the environmental time(perhaps due to rapid travel across multiple time zones) the optimal control input is computed forvarious choices of ∆ t with the adaptive reduction framework (from Equations (14), (19) and (20))and the first order accurate phase-amplitude reduction strategy (from Equations (21) and (24)).For when using the first order phase-amplitude reduction strategy, β is set to zero since the nominal24-hour light-dark cycle makes it difficult to limit the isostable coordinates. For many choices of ∆ t ,the phase-only optimal control equations (27) have solutions that grow to infinity in finite time;consequently this optimal control formulation will not be considered for this application. In alloptimal control computations, t = 0 which corresponds to the middle of subjective night. Whenpossible, T f is taken to be 24 hours. For positive time shifts and particularly large magnitudenegative time shifts, solutions of the optimal control equations were not able to be found using T f = 24 hours. In these cases, T f is instead taken to be 48 hours.The numerically determined optimal control inputs are applied to the full model equations (30)with results shown in Figure 7. The recovery time in panel A is defined to be the time required(after the time shift by ∆ t of the light-dark cycle) for the phase to return to and stay within twohours of its steady state behavior. The nominal recovery time that results when ∆ L is zero is alsoshown for reference. Note that the moment that F ( t ) crosses 0.045 is used to correspond to thetime that θ = 0 and that this can only be observed once per cycle – for this reason, if the phase is18igure 7: Panel A shows recovery times when applying inputs obtained from the adaptive reductionand first order accurate phase-isostable reduction are shown as blue and red lines, respectively. Forreference, the black line show the uncontrolled recovery times. Recovery is substantially hastened foroptimal inputs obtained from the adaptive reduction. Conversely, it is often worse to apply inputsobtained from the first order accurate phase-isostable reduction than it is to do nothing. PanelsB and C show example optimal control inputs for a positive and negative time shift, respectively.The stimuli in panel B are qualitatively similar for each optimal control both suggesting that twoshort pulses of light will engender the required + 7 hour phase shifts. Nevertheless, the differencesin timing and magnitude of these stimuli yield qualitatively different behaviors when applied to thefull model (30) as shown in panel D. Similarly, the stimuli shown in panel C are nearly identical for t ∈ [0 , θ = 0 after the optimal stimulus is applied, the recovery time istaken to be T f . For each optimal control determined from the adaptive reduction, recovery occursat or very close to to T f . The optimal control derived from the first order accurate phase-isostablereduction speeds recovery time for phase advances. For phase delays, however, recovery is actuallydelayed when the optimal control identified from the first order phase-isostable reduction is used.Panel B (resp. C) shows control inputs obtained for a ∆ t = +7 hour (resp. -10 hour) time shift. Theresulting model outputs when stimulus is applied to the full model equations (30) are shown in panelD (resp. E). In both cases, the input derived from the adaptive reduction is qualitatively similarto the input from the input derived using the phase-isostable reduction, however, the quantatativedifferences between each input result in substantially different recovery times.Finally, the individual cell dynamics during recovery are considered in Figure 8. For a timeshift ∆ t = +10 hours that occurs instantaneously starting at t = 0 hours, traces of E i ( t ) for30 representative oscillators are shown as colored lines in addition to F ( t ) (solid black line) anda reference F ( t ) that is fully entrained to the +10 hour time shifted light-dark cycle. Each dotof the raster plot in Panel B represents the moment that a single cell in the simulation crossesthe threshold E i ( t ) = 0 . L ( t ) is only applied during t ∈ [0 , r rot is the19igure 8: Panel A shows traces of E i ( t ) for 30 representative oscillators from the larger populationwhen a time shift of ∆ t = +10 hours is applied starting at t = 0. The optimal control (panel C) isobtained using the adaptive reduction and is applied in order to rapidly reentrain the populationoscillation to the new 24-hour light-dark cycle. Each dot from the raster plot in panel B representsthe time that E i crosses the threshold 0.035 with a positive slope taken to correspond to θ = 0 for theindividual oscillators. Panel D shows the rotating order parameter from (32) during the simulation.Panels E-H provide analogous information for a time shift of ∆ t = − r rot takes paths through the center in the complex plane.Kuramoto order parameter [27] taken in a rotating reference frame: r rot ( t ) = 1 N N (cid:88) k =1 e i ( θ k ( t ) − πt/ , (32)where θ i ( t ) is the phase of oscillator i at time t . Here, r rot will be referred to as the rotating orderparameter. Intuitively, Equation (32) contains similar information to the standard definition of theorder parameter; | r rot | gives a sense of the cohesion of the oscillators, where | r rot | = 1 when thenetwork is perfectly synchronized and | r rot | = 0 when the neurons are in a splay state. The term2 πt/
24 is subtracted from the phase in (32) so that when a time shift of ∆ t is applied to the model,the argument of the complex number r rot ( t ) shifts by 2 π ∆ t/
24 in steady state (i.e., once the modelis fully reentrained to the shifted 24-hour light-dark cycle). Because information about the phasesis not directly available, it is assumed that each oscillator reaches a phase of θ = 0 when its valueof E i crosses 0.035 with a positive slope, and all other phase values are linearly interpolated fortimes in between. Panels E-H show similar information for a time shift of ∆ t = − Discussion and Conclusion
In this work, an optimal control strategy is developed for manipulating oscillation timing in thepresence of a large magnitude entraining stimulus. In contrast to other recently proposed strategiesthat leverage phase reduced models for control of oscillation timing [34], [37], [35], [58] the frame-work considered here is applicable for arbitrary large magnitude inputs. Consequently, larger shiftsin phase can be considered without sacrificing accuracy. While this work is primarily motivated bythe search for control protocols that hasten recovery from jet-lag in order to mitigate its negativeside effects, the control formulation is quite general and could be implemented in other applicationswhere rapid modification of oscillation timing is required.In order to accurately account for the influence of large magnitude inputs in a reduced or-der setting, the recently developed phase-amplitude reduction framework is employed here. Byconsidering a family of limit cycles that emerge for different parameter sets, the nominal systemparameters can be adaptively updated in order to limit the deviation from the nominal limit cycle,thereby minimizing error in the reduced order equations. In general, the dimension of the adaptivereduction grows with both the number of isostable coordinates required and the number of theadaptive parameters. However, as shown in this work, provided the amplitude corrections can becaptured by a single isostable coordinate and Q ( θ, p ) is nonzero for all allowable θ and p , the adap-tive parameter can be updated in such a way that the isostable coordinate dynamics are eliminated,yielding a reduced order set of equations with only two dimensions. A subsequent optimal controlformulation for manipulating the timing of oscillations based on Pontryagin’s minimum princi-ple [24] is considered in order to make comparisons with previously published results, however, awide variety of optimal control formulations would be possible given the resulting two-dimensionalreduced order representation of the dynamics.For both the simple three-dimensional oscillator model (28) and the model comprised of 3000heterogeneous, noisy coupled oscillators (30) the optimal control formulation using the adaptivereduction strategy vastly outperforms a strategy that uses a phase-only reduction [34] and a strategythat uses a first order accurate phase-amplitude reduction [35]. When no entraining stimulus isconsidered (as is the case in results presented in Figure 2), either reduction framework is sufficientfor identifying stimuli that can shift the oscillation timing. However, the phase-only and first orderphase-amplitude reduction strategy begins to degrade rapidly as the magnitude of the mandatedshift grows. Such degradation is not observed when the adaptive reduction is used, even whenconsidering the largest possible phase shifts of ∆ t = ±
12 hours (i.e., ∆ θ ≈ ± π ). The incorporationof a large magnitude input that simulates the influence of a nominal 24-hour external light-darkcycle in the optimal control formulation renders the results when phase-only and first order accuratephase-amplitude reduced equations unusable; resulting recovery times in Figure 3 are similar to theuncontroled recovery for these formulations. Conversely, optimal control waveforms obtained whenusing the adaptive reduction perform as desired with negligible errors. A similar story emerges inthe results from Figure 7 for the coupled population model (30).A secondary focus of this work considers the identification of the necessary terms of an adaptivephase reduction using data-driven techniques. The direct method [21], [38] is a well-established,data-driven technique for inferring phase reduced equations of the form (1) when the underlyingdynamical equations are unknown (e.g., in experimental systems). While related methods have beenproposed to infer numerical equations that characterize the behavior of the amplitude coordinates[60], [56], there is still much work to be done. As illustrated for the coupled population of limitcycle oscillators (30), provided that the adaptive system parameter is chosen appropriately, all ofthe terms of the adaptive phase reduction can be determined by applying a series of pulse and stepfunction inputs and inferring the resulting change in the phase and isostable coordinates. While21his process requires more data than simply applying the direct method to compute the phasereduced equations (1), ignoring the extra information afforded by the adaptive reduction (e.g., byonly implementing the first order accurate phase-amplitude reduction in the results shown in Figure7) renders the resulting optimal control outputs useless.When considering the problem of recovery from circadian misalignment, an interesting resultobserved here is that for both large magnitude time advances and time delays, the optimal controlwaveforms first desynchronize the individual oscillators within the population, effectively destroy-ing the collective oscillation, before resynchronizing with the target phase of oscillation. Such anotion is related to the observation from in vito experiments that the cells of the suprachiasmaticnucleus reentrain faster to shifts in environmental time if they are first desynchronized prior to thetime shift [2]. While the level of desynchronization was not explicitly considered by the optimalcontrol formulation in Section 5, the resulting optimal control inputs significantly desynchronizethe population during the phase shift. In future work, it would be of interest to investigate themathematical mechanisms by which desynchronization allows for more rapid shifts in the popu-lation phase of oscillation and whether these observations could be leveraged to develop bettertreatments for rapid recovery from circadian misalignment.While the adaptive phase reduction used here shows promise for control of oscillation timing inapplications where large magnitude inputs are required, there are many opportunities for improve-ment. From a theoretical perspective, little work has been done concerning the design of function G p of the adaptive reduction (11) that is required to limit the magnitude of isostable coordinates ψ , . . . , ψ β . Of course, when only a small number of isostable coordinates are considered, as is thecase in this work, adequate choices for G p are relatively easy to identify. The same is not true whenmany isostable coordinates are considered and more consideration about the design of G p whenmultiple isostable coordinates are required would be warranted. As a matter of practical imple-mentation, it would be useful to identify strategies to experimentally infer the necessary terms ofthe adaptive reduction from (11) that allow for arbitrary adaptive parameters to be considered. Asignificant difficulty in inferring the terms of the adaptive reduction using data-driven methods is inthe determination of U e ( t, p, x ). Note that the definition from Equation (7), U e requires knowledgeof the underlying dynamical equations which are usually unknown in experimental applications.The explicit assumption used in this work that U ( t ) = δu ( t ), and F ( x, p ) = F ( x ) + pδ where δ ∈ R N (i.e., that a shift in the nominal parameter p provides the same effect as a shift in u ( t ))provides a workaround to this limitation that allows U e to be computed exactly. Nevertheless, thisassumption severely limits the choice of the adaptive parameters and it would be useful to identifystrategies that can infer the terms of the adaptive reduction for more general choices of the adaptiveparameter. Finally, the amount of data required to infer the terms of the adaptive reduction usingthe methodology suggested in Section 5.1 may be particularly onerous in some applications – eachparameter set considered requires a sufficient number of trials performed using both pulse and stepfunction inputs to provide adequate data to fit functions to z , i , D , and Q as illustrated in Figure5. Strategies that can be used to infer the terms of the adaptive reduction more efficiently, suchas those that use deep learning approaches [31], [65] would certainly be of interest and will beinvestigated in future work.This material is based upon the work supported by the National Science Foundation (NSF)under Grant No. CMMI-1933583. 22 eferences [1] D. M. Abrams and S. H. Strogatz. Chimera states for coupled oscillators. Physical ReviewLetters , 93(17):174102, 2004.[2] S. An, R. Harang, K. Meeker, D. Granados-Fuentes, C. A. Tsai, C. Mazuski, J. Kim, F. J.Doyle, L. R. Petzold, and E. D. Herzog. A neuropeptide speeds circadian entrainment by reduc-ing intercellular synchrony.
Proceedings of the National Academy of Sciences , 110(46):E4355–E4361, 2013.[3] N. Bagheri, J. Stelling, and F. J. Doyle III. Circadian phase resetting via single and multiplecontrol targets.
PLoS Computational Biology , 4(7):e1000104, 2008.[4] P. Benner, S. Gugercin, and K. Willcox. A survey of projection-based model reduction methodsfor parametric dynamical systems.
SIAM Review , 57(4):483–531, 2015.[5] P. C. Bressloff and J. N. MacLaurin. A variational method for analyzing stochastic limit cycleoscillators.
SIAM Journal on Applied Dynamical Systems , 17(3):2205–2233, 2018.[6] E. Brown, P. Holmes, and J. Moehlis. Globally coupled oscillator networks. In
Perspectivesand Problems in Nonlinear Science , pages 183–215. Springer, 2003.[7] M. Budiˇsi´c, R. Mohr, and I. Mezi´c. Applied Koopmanism.
Chaos: An Interdisciplinary Journalof Nonlinear Science , 22(4):047510, 2012.[8] O. Castej´on, A. Guillamon, and G. Huguet. Phase-amplitude response functions for transient-state stimuli.
J. Math. Neurosci , 3:13, 2013.[9] J. Cui, C. C. Canavier, and R. J. Butera. Functional phase response curves: A method forunderstanding synchronization of adapting neurons.
Journal of Neurophysiology , 102(1):387–398, 2009.[10] D. A. Dean, D. B. Forger, and E. B. Klerman. Taking the lag out of jet lag through model-basedschedule design.
PLoS Computational Biology , 5(6):e1000418, 2009.[11] C. O. Diekman and A. Bose. Entrainment maps: a new tool for understanding properties ofcircadian oscillator models.
Journal of Biological Rhythms , 31(6):598–616, 2016.[12] C. O. Diekman and A. Bose. Reentrainment of the circadian pacemaker during jet lag: east-west asymmetry and the effects of north-south travel.
Journal of Theoretical Biology , 437:261–285, 2018.[13] F. D¨orfler, M. Chertkov, and F. Bullo. Synchronization in complex oscillator networks andsmart grids.
Proceedings of the National Academy of Sciences , 110(6):2005–2010, 2013.[14] G. B. Ermentrout and D. H. Terman.
Mathematical Foundations of Neuroscience , volume 35.Springer, New York, 2010.[15] Gr´egory G. Dumont, G. B. Ermentrout, and B. Gutkin. Macroscopic phase-resetting curvesfor spiking neural networks.
Physical Review E , 96(4):042311, 2017.[16] E. Genge, E. Teichmann, M. Rosenblum, and A. Pikovsky. High-order phase reduction forcoupled oscillators. arXiv preprint 2007.14077 , 2020.2317] D. Gonze, S. Bernard, C. Waltermann, A. Kramer, and H. Herzel. Spontaneous synchronizationof coupled circadian oscillators.
Biophysical Journal , 89(1):120–129, 2005.[18] J. Guckenheimer. Isochrons and phaseless sets.
Journal of Mathematical Biology , 1(3):259–273,1975.[19] P. Holmes, J. L. Lumley, G. Berkooz, and C. W. Rowley.
Turbulence, coherent structures,dynamical systems and symmetry . Cambridge University Press, New York, 1996.[20] A. B. Holt, D. Wilson, M. Shinn, J. Moehlis, and T. I. Netoff. Phasic burst stimulation: aclosed-loop approach to tuning deep brain stimulation parameters for parkinson’s disease.
PLoS Computational Biology , 12(7):e1005011, 2016.[21] E. M. Izhikevich.
Dynamical Systems in Neuroscience: The Geometry of Excitability andBursting . MIT Press, London, 2007.[22] E. Kaiser, N. J. Kutz, and S. L. Brunton. Data-driven discovery of Koopman eigenfunctionsfor control. arXiv preprint arXiv:1707.01146 , 2017.[23] Y. Kawamura, H. Nakao, K. Arai, H. Kori, and Y. Kuramoto. Collective phase sensitivity.
Physical Review Letters , 101(2):024101, 2008.[24] D. Kirk.
Optimal Control Theory . Dover Publications, New York, 1998.[25] T. W. Ko and G. B. Ermentrout. Phase-response curves of coupled oscillators.
Physical ReviewE , 79(1):016211, 2009.[26] N. J. Kopell, H. J. Gritton, M. A. Whittington, and M. A. Kramer. Beyond the connectome:the dynome.
Neuron , 83(6):1319–1328, 2014.[27] Y. Kuramoto.
Chemical Oscillations, Waves, and Turbulence . Springer-Verlag, Berlin, 1984.[28] W. Kurebayashi, S. Shirasaka, and H. Nakao. Phase reduction method for strongly perturbedlimit cycle oscillators.
Physical Review Letters , 111(21):214101, 2013.[29] B. Letson and J. E. Rubin. LOR for analysis of periodic dynamics: A one-stop shop approach.
SIAM Journal on Applied Dynamical Systems , 19(1):58–84, 2020.[30] Z. Levnaji´c and A. Pikovsky. Phase resetting of collective rhythm in ensembles of oscillators.
Physical Review E , 82(5):056202, 2010.[31] B. Lusch, N. J. Kutz, and S. L. Brunton. Deep learning for universal linear embeddings ofnonlinear dynamics.
Nature Communications , 9(1):1–10, 2018.[32] A. Mauroy, I. Mezi´c, and J. Moehlis. Isostables, isochrons, and Koopman spectrum for theaction–angle representation of stable fixed point dynamics.
Physica D: Nonlinear Phenomena ,261:19–30, 2013.[33] I. Mezi´c. Spectrum of the Koopman operator, spectral expansions in functional spaces, andstate-space geometry.
Journal of Nonlinear Science , pages 1–55, 2019.[34] J. Moehlis, E. Shea-Brown, and H. Rabitz. Optimal inputs for phase models of spiking neurons.
ASME Journal of Computational and Nonlinear Dynamics , 1(4):358–367, 2006.2435] B. Monga and J. Moehlis. Optimal phase control of biological oscillators using augmentedphase reduction.
Biological Cybernetics , 113(1-2):161–178, 2019.[36] R. Y. Moore, J. C. Speh, and R. K. Leak. Suprachiasmatic nucleus organization.
Cell andTissue Research , 309(1):89–98, 2002.[37] A. Nabi, T. Stigen, J. Moehlis, and T. Netoff. Minimum energy control for in vitro neurons.
Journal of Neural Engineering , 10(3):036005, 2013.[38] T. Netoff, M. A. Schwemmer, and T. J. Lewis. Experimentally estimating phase response curvesof neurons: theoretical and practical issues. In
Phase Response Curves in Neuroscience , pages95–129. Springer, 2012.[39] E. Ott and T. M. Antonsen. Low dimensional behavior of large systems of globally coupledoscillators.
Chaos: An Interdisciplinary Journal of Nonlinear Science , 18(3):037113, 2008.[40] Y. Park and B. Ermentrout. Weakly coupled oscillators in a slowly varying world.
Journal ofComputational Neuroscience , 40(3):269–281, 2016.[41] B. Pietras and A. Daffertshofer. Network dynamics of coupled oscillators and phase reductiontechniques.
Physics Reports , 2019.[42] K. Pyragas and V. Noviˇcenko. Phase reduction of a limit cycle oscillator perturbed by a strongamplitude-modulated high-frequency force.
Physical Review E , 92(1):012910, 2015.[43] S. M. Reppert and D. R. Weaver. Coordination of circadian timing in mammals.
Nature ,418(6901):935, 2002.[44] M. Rosenblum and A. Pikovsky. Numerical phase reduction beyond the first order approxi-mation.
Chaos: An Interdisciplinary Journal of Nonlinear Science , 29(1):011105, 2019.[45] F. Della Rossa, L. Pecora, K. Blaha, A. Shirin, I. Klickstein, and F. Sorrentino. Symmetriesand cluster synchronization in multilayer networks.
Nature Communications , 11(1):1–17, 2020.[46] R. L. Sack. Jet lag.
New England Journal of Medicine , 362(5):440–447, 2010.[47] J. T. C. Schwabedal and A. Pikovsky. Phase description of stochastic oscillations.
PhysicalReview Letters , 110(20):204102, 2013.[48] K. Serkh and D. B. Forger. Optimal schedules of light exposure for rapidly correcting circadianmisalignment.
PLoS Computational Biology , 10(4):e1003523, 2014.[49] S. Shirasaka, W. Kurebayashi, and H. Nakao. Phase-amplitude reduction of transient dynamicsfar from attractors for limit-cycling systems.
Chaos: An Interdisciplinary Journal of NonlinearScience , 27(2):023119, 2017.[50] S. H. Strogatz, D. M. Abrams, A. McRobie, B. Eckhardt, and E. Ott. Theoretical mechanics:Crowd synchrony on the millennium bridge.
Nature , 438(7064):43, 2005.[51] K. Taira, S. L. Brunton, S. T. M. Dawson, C. W. Rowley, T. Colonius, B. J. McKeon, O. T.Schmidt, S. Gordeyev, V. Theofilis, and L. S. Ukeiley. Modal analysis of fluid flows: anoverview.
AIAA Journal , pages 4013–4041, 2017.2552] P. A. Vela.
Averaging and control of nonlinear systems . PhD thesis, California Institute ofTechnology, 2003.[53] K. C. A. Wedgwood, K. K. Lin, R. Thul, and S. Coombes. Phase-amplitude descriptions ofneural oscillator models.
The Journal of Mathematical Neuroscience , 3(1):2, 2013.[54] D. Wilson. An adaptive phase-amplitude reduction framework without O ( (cid:15) ) constraints oninputs. Submitted .[55] D. Wilson. Isostable reduction of oscillators with piecewise smooth dynamics and complexFloquet multipliers.
Physical Review E , 99(2):022210, 2019.[56] D. Wilson. A data-driven phase and isostable reduced modeling framework for oscillatorydynamical systems.
Chaos: An Interdisciplinary Journal of Nonlinear Science , 30(1):013121,2020.[57] D. Wilson. Phase-amplitude reduction far beyond the weakly perturbed paradigm.
PhysicalReview E , 101(2):022220, 2020.[58] D. Wilson and B. Ermentrout. Greater accuracy and broadened applicability of phase reductionusing isostable coordinates.
Journal of Mathematical Biology , 76(1-2):37–66, 2018.[59] D. Wilson and B. Ermentrout. An operational definition of phase characterizes the transientresponse of perturbed limit cycle oscillators.
SIAM Journal on Applied Dynamical Systems ,17(4):2516–2543, 2018.[60] D. Wilson and B. Ermentrout. Augmented phase reduction of (not so) weakly perturbedcoupled oscillators.
SIAM Review , 61(2):277–315, 2019.[61] D. Wilson and B. Ermentrout. Phase models beyond weak coupling.
Physical Review Letters ,123(16):164101, 2019.[62] D. Wilson, S. Faramarzi, J. Moehlis, M. R. Tinsley, and K. Showalter. Synchronization ofheterogeneous oscillator populations in response to weak and strong coupling.
Chaos: AnInterdisciplinary Journal of Nonlinear Science , 28(12):123114, 2018.[63] D. Wilson and J. Moehlis. Isostable reduction of periodic orbits.
Physical Review E ,94(5):052213, 2016.[64] A. Winfree.
The Geometry of Biological Time . Springer Verlag, New York, second edition,2001.[65] E. Yeung, S. Kundu, and N. Hodas. Learning deep neural network representations for koopmanoperators of nonlinear dynamical systems. In , pages 4832–4839. IEEE, 2019.[66] A. Zlotnik, R. Nagao, I. Z. Kiss, and J. S. Li. Phase-selective entrainment of nonlinear oscillatorensembles.