PHASE SENSITIVITY ANALYSIS OF LIMIT CYCLE OSCILLATORS

A. K. Wilkins1, B. Tidor23, J. White3, P. I. Barton1*
1Dept. of Chemical Engineering
2Biological Engineering Division
3Dept. of Electrical Engineering and Computer Science
Massachusetts Institute of Technology, USA
*pib@mit.edu

A comprehensive analysis is presented of the different types of phase sensitivities of limit cycle oscillators with respect to parameters. A general framework based on the solution of a boundary value problem is presented which permits exact calculation of the desired quantities. Furthermore, a decomposition of the state sensitivities into a sum of three terms each with a specific contribution to the overall sensitivity is discussed in detail. Numerical examples are used to demonstrate the techniques and validate the results. Example systems are taken from the biological literature, where limit cycle oscillators of different levels of detail are used to model the circadian clock mechanism. 

The sensitivity analysis of limit cycle oscillators is, compared to the sensitivity analysis of other dynamical systems, complicated by several factors. First, from any initial condition within the region of attraction, the periodic orbit is reached only in asymptotic fashion, but it is of interest to study its exact properties. Furthermore, it is desired to analyze derived properties which are implicit in the mathematical description of the limit cycle, such as amplitude, period and phase of oscillation, and their parametric dependencies.
A boundary value formulation, similar to the one proposed by Rosenwasser and Yusupov [1], avoids dealing with transients and allows the direct computation of the period sensitivities of the limit cycle. A boundary value problem (BVP) is formulated similarly to those used in circuit analysis [2], in such a fashion that the first n equations ensure the periodicity of the solution, and the (n+1)st equation is a phase locking condition (PLC) which defines a reference for time zero. When these equations are formally differentiated with respect to the parameters, the resulting matrix equation can be solved for the period sensitivities and the sensitivities of the initial conditions on the limit cycle, as found from the solution of the BVP. Using these initial condition sensitivities, the sensitivity trajectories can be integrated along with the dynamic system, and amplitude sensitivities can be obtained from the sensitivities at the extrema of the various state variables [1].
The work presented here shows that the choice of PLC is of significance not only for the numerical well-posedness of the resulting BVP, but more importantly for the computation and interpretation of meaningful sensitivities of derived properties, in particular those of the peak-to-peak sensitivities, as described below. 
As shown in [1,3-4], the state sensitivities can be separated into an unbounded, period-dependent part and a strictly periodic part,

					(1)
where  are the sensitivities of the state variables with respect to the parameters, are the period sensitivities with respect to the parameters, and are the -periodic sensitivities of the state variables with respect to the parameters at a fixed period. 
This work introduces an exact and unique decomposition of the periodic part  into a periodic, phase dependent part and a periodic, amplitude dependent part. (A similar decomposition into amplitude and phase components was used for oscillator noise analysis [5]). The full sensitivity is thus composed of three specific contributions:

whereis a constant row vector that describes the influence of the PLC on the solution of the BVP and contains information about the sensitivities in the direction tangential to the periodic orbit. are the state sensitivities for constant period and constant phase, in other words theterm describes only the influence of the parameters on the amplitudes of the different state variables and therefore the shape of the limit cycle.

This new decomposition is the starting point for the definition and computation of different types of phase sensitivities. First, a sensitivity of the phase set by the PLC with respect to the parameters is discussed. This quantity is then used to define a particular type of sensitivity, which is unique to oscillatory systems, and is termed a “peak-to-peak” sensitivity. It describes how an infinitesimal change in parameters would cause a change in the time between two peaks or troughs of different state variables. The choice of PLC is obvious in this case, an extremum of a state variable of interest is chosen as the time reference. A simple equation is derived for the exact computation of these sensitivities and the results are compared to the proportional sensitivity of the peak-to-peak distance with respect to overall period change (i.e., stretching or compressing of the oscillation). In doing so, information about the flexibility of the system is gained, in the sense of whether the time span between specific events (e.g., peaks) can be varied independently of the overall period. 
This quantity is of interest in the context of circadian biology; while the total length of day is constant at 24 hours, the length of the sunlit day, or the time between dawn and dusk, undergoes seasonal variation. It could be hypothesized that the clock mechanism might have the capability of adjusting to this variation by varying the peak-to-peak distances of pertinent state variables without varying the period of oscillation. This hypothesis is evaluated using peak-to-peak sensitivity analysis. 

A different type of phase sensitivity analysis is the analysis of phase shifts caused by perturbing the state variables from the limit cycle. Such perturbations cause temporary deviations from the periodic orbit that decay over time, but result in permanent phase shifts that are a functions of the perturbation and the system parameters. The sensitivity of this permanent phase shift with respect to the perturbed variables and also with respect to the system parameters is derived and computed exactly and efficiently. 
The sensitivity of the phase shift with respect to the perturbed states was approximated in the past using a limit finding process [4]. We show, however, that this sensitivity information is contained in the first left eigenvector of the Monodromy matrix of the limit cycle oscillator, with appropriate scaling. 

The sensitivities described above were computed for two model systems. The first test system is a simple model of the circadian clock, called the Goodwin oscillator, composed of 3 state variables and 6 parameters [6]. The second test system is the most detailed ordinary differential equation model of the circadian clock currently published, which consists of 73 state variables and 231 parameters [7]. All results were compared to finite difference approximations, to demonstrate the methods are accurate yet far more efficient. 

References

[1] Rosenwasser E, Yusupov R: Sensitivity of Automatic Control Systems, CRC Press, Boca Raton, FL 2000.

[2] Kundert KS, White JK, Sangiovanni-Vincentelli AL: Steady-State Methods for Simulating Analog and Microwave Circuits, Kluwer Academic Publishers, Boston, MA 1990. 

[3] Larter R, Rabitz H, Kramer M: Sensitivity analysis of limit cycle oscillators with application to the Brusselator. J Chem. Phys 80, 4120 – 4218 (1984)

[4] Kramer MA, Rabitz H, Calo JM: Sensitivity Analysis of oscillatory systems. Appl Math Modelling 8, 328-340 (1984)

[5] Kaertner F X: Analysis of white and f - noise in oscillators. Int J circuit Theory Appl , 485-519 (1990)

[6] Ruoff P, Rensing L: The Temperature-Compensated Goodwin Model Simulates Many Circadian Clock Properties J theor Biol 179, 275 – 285 (1996)

[7] Forger D B and Peskin CS: A detailed predictive model of the mammalian circadian clock. Proc Natl Acad Sci U S A 100, 14806 - 14811 (2003)