Pacific Symposium on Biocomputing 13:402-413(2008) September 26, 2007 22:2 WSPC - Proceedings Trim Size: 9in x 6in agarwal E FFI C I E N T M U LT I S C A L E S I M U L AT I O N O F C I R C A D I A N R H Y T H M S U S I N G AU T O M AT E D P H A S E M AC RO M O D E L L I N G T E C H N I Q U E S S H ATA M AG A RWA L Indian Institute of Technology, Kanpur shatam@iitk.ac.in JA I J E E T ROY C H OW D H U RY University of Minnesota, Twin Cities jr@umn.edu Circadian rhythm mechanisms involve multi-scale interactions between endogenous DNA-transcription oscillators. We present the application of efficient, numerically wellconditioned algorithms for abstracting (potentially large) systems of differential equation models of circadian oscillators into compact, accurate phase-only macromodels. We apply and validate our auto-extracted phase macromodelling technique on mammalian and Drosophila circadian systems, obtaining speedups of 9 - 13× over conventional timecourse simulation, with insignificant loss of accuracy, for single oscillators being synchronized by day/night light variations. Further, we apply the macromodels to simulate a system of 400 coupled circadian oscillators, achieving speedups of 240× and accurately reproducing synchronization and locking phenomena amongst the oscillators. We also present the use of parameterized phase macromodels for these circadian systems, and elucidate insights into circadian timing effects directly provided by our auto-extracted macromodels. 1. Introduction Circadian rhythms are amongst the most fundamental of physiological processes. They are found in virtually all organisms, ranging from unicellular (e.g., amoebę, bacteria) to complex multicellular higher organisms (e.g., human beings). These daily rhythms, of period about 24 hours, are associated with periodic changes in hormones controlling sleep/wakefulness, body temperature, blood pressure, heart rate and other physiological variables. Importantly, circadian rhythms are endogenous or autonomous; however, they are typically influenced by external cues, such as light. Progress in quantitative biology has established that such rhythms stem fundamentally from the molecular level,1,2 involving complex chains of biochemical reactions featuring a number of key proteins/hormones (such as melatonin and melanopsin), whose levels rise and fall during the course of the day. These biochemical reactions, which take place both within individual cells and at an extracellular level, function as biological oscillators or body clocks.3 Quantitative understanding, simulation and control of circadian rhythms is of great practical importance. Applications include devising medical remedies for rhythm disorders (e.g., insomnia, fatigue, jet lag, etc.), synthetic biology (where a goal is to "program" artificial rhythms that are biologically viable), artificially extending periods of wakefulness/alertness (e.g., for military purposes), and so on. Improved understanding of circadian rhythm mechanisms has led to increased awareness of how pervasively they affect virtually every aspect of the life of an organism. Hence, their simulation/analysis is an important endeavour in the biological domain.1,2 Although individual oscillators constitute the fundamental core of circadian rhythm mechanisms, the rich circadian functionality of multicellular organisms results from the interactions of many oscillators over multiple temporal and spatial scales. Observations of periodicity in behavior, metabolism, body temperature, etc., indicate that coupling/coherence mechanisms play a key role. Hierarchical organization of the circadian system, from the fundamental DNA transcription/translation level to endocrine system levels, involves complex oscillator interactions. The complex connectivity and high dimensionality of such coupled oscillator networks, which lead to unique effects such as synchronization and injection locking/pulling,4,5 make them difficult to understand at the intuitive or analytical level, thus engendering the need for efficient and powerful simulation and analysis tools with multiscale capabilities. Several oscillatory mathematical models are available for circadian rhythms1,2 that capture the dynamics of the relevant molecular biochemical reactions (see Section 2 for details). These models are in the form of systems of differential-algebraic equations (DAEs) or ordinary differential equations (ODEs). The prevalent technique today for their simula- Pacific Symposium on Biocomputing 13:402-413(2008) September 26, 2007 22:2 WSPC - Proceedings Trim Size: 9in x 6in agarwal tion is to run initial value simulations. While such "time-course integration" of ODEs/DAEs has the advantage of generality, it suffers from serious disadvantages for oscillators, which are inherently marginally stable.6 For initial-value simulations, marginally stable systems tend to require orders of magnitude more computation for a specified accuracy, particularly phase/timing accuracy, than stable systems; even for individual oscillators, very small timesteps (e.g., many hundreds per oscillation cycle) are typically needed, leading to high computational cost. The situation worsens for coupled oscillator systems, which typically feature multiple time scales; e.g., envelopes7,8 typically feature much longer time scales than individual oscillation cycles. In electronic circuit design, automated nonlinear phase macromodel extraction techniques5,6,9 have proven effective in solving sucsuch oscillatory problems. Given any oscillator as a system of DAEs or ODEs (however complicated), efficient and well-conditioned numerical techniques extract a scalar nonlinear differential equation, the phase macromodel. This macromodel captures the dynamic response of the oscillator's phase or timing characteristics to external influences. It has been shown that such "PPV" (Perturbation Projection Vector) phase macromodels are able to accurately capture the gamut of phase/frequency-related dynamics of oscillators; most importantly locking, synchronization and phase noise (timing jitter) effects.6,10 Using the PPV macromodel instead of the original DAEs/ODEs confers important advantages: large simulation speedups due to system size reduction, the ability to use larger timesteps than for the original system, abstraction to the phase or timing level, precise insight about timing influences without the need for simulation, etc.. These advantages are especially pronounced for systems of many coupled oscillators spanning different temporal and spatial scales.11 In this work, we present the first application of PPV-based automated nonlinear timeshifted macromodelling methods to biological systems, focussing on circadian rhythms. We use PPV phase macromodels6,12 to model circadian oscillators and show that they are considerably more efficient than standard "time course" simulations. PPV models alleviate the lack of accuracy and general applicability of a widely used prior phase model (Kuramoto's model, see below), while retaining its advantage of relative simplicity and computational efficiency. PPVs provide direct insights into the effects of external stimuli, such as slowing down/speeding up of circadian rhythms; for example, it is easy to determine when and how to apply a light pulse for greatest de-synchronization. Using PPV macromodels, we are able to efficiently produce plots of circadian lock range vs amplitude of external stimuli; this is valuable for guiding experiments, explaining observations, and designing new ("synthetic") DNA/protein based biological clock networks. Furthermore, we also present the application of parameterized PPV macromodels,13 which directly incorporate effects of parameter changes into our phase-only models of circadian rhythms. Being able to directly and quantitatively predict the impact of parameter changes on phase, frequency and timing behaviour is of significant biological value. Indeed, the PPV constitutes a rigorous, Floquet-theoretic generalization of Winfree's seminal concept of timing maps and phase sensitivity functions,14 used within "phase-only models" of oscillators popular in computational biology. Kuramoto15 applied the theory of asymptotics to find sinusoidal expressions for Winfree's phase sensitivity functions; these are widely used for phase-only models, since they are capable of capturing synchronization effects. However, the sinusoidal simplifications inherent in Kuramoto's model lead to significant inaccuracies for non-sinusoidal oscillators (including circadian ones). These inaccuracies compound in large networks and often lead to qualitatively incorrect conclusions about, e.g., collective synchronization.11,16 From a utilitarian point of view, the usefulness of PPV macromodels over Kuramoto is twofold: firstly, PPVs represent the correct (often highly non-sinusoidal) phase sensitivity functions of the oscillator; secondly, the PPV macromodel is generated via an algorithmic computational procedure from the oscillator's DAE/ODE description, typically taking seconds or less to reduce systems 1000s of equations in size. Pacific Symposium on Biocomputing 13:402-413(2008) September 26, 2007 22:2 WSPC - Proceedings Trim Size: 9in x 6in agarwal We apply PPV macromodels to two different circadian rhythm models: one for mammals1,17 and one for Drosophila melanogaster (the fruit fly, shown in Fig. 1).2 We show that the PPV macromodels are significantly faster to simulate than the original equation systems even for single oscillators (9x speedup for the mammalian clock, and 13x speedup for the Drosophila clock). Modelling light as an external input impinging upon circadian oscillators, we con- Fig. 1: Male Drosophila firm injection locking using PPV macromodels and obtain plots (fruit fly). of lock range vs amplitude of the light signal. We comment on the biological significance of the shapes of important components of the PPV. We then use PPV macromodels to rapidly explore synchronization behaviour in a network of 400 coupled oscillators, obtaining speedups of about 240× over standard time-course simulation. Finally, using parameterized PPV macromodels, we predict the effects of varying a number of model parameters on oscillation frequency and lock range. The remainder of the paper is organized as follows. In Section 2.1, we provide background on circadian rhythms and their mechanisms, followed by a review of mathematical models for circadian rhythms in Section 2.2. Oscillator and phase macromodels are then introduced in Section 3; a brief review of PPV macromodels, injection locking analysis and parameterized PPVs is provided. Finally, in Section 4, results and speedups are presented. 2. Background and Previous Work 2.1. Circadian rhythms Circadian rhythms are generated by "clock genes", which encode genetic instructions that produce certain proteins whose levels oscillate during the course of the day. These oscillating biochemical signals control various functions, such as sleep/waking cycles ­ in other words, they constitute our "internal biological clock", which adapts to the daily cycle of day and night. However, the natural period of this internal clock is not exactly 24 hours; it is typically longer if the organism is kept isolated and away from external cues,18 most importantly light (these cues are called Zeitgebers). Therefore, the internal clock needs to be "reset" every day, in order to keep the organism's bodily rhythms synchronized with the external world's day/night cycle. Higher organisms are often composed of billions of cells. The nucleus of each cell contains the genetic material DNA, a long chain-like linear molecule built up of many links. RNA, also a nucleic acid polymer, serves as a DNA template for the translation of genes into proteins. The process of formation of an RNA molecule from a particular DNA is called transcription. Unlike DNA, RNA is capable of leaving the nucleus and moving into the cytoplasm. There, with the help of enzymes, specific RNA strings get converted to specific proteins responsible for different bodily functions. Some of these proteins return to the nucleus, forming complexes by binding to other proteins, some of which inhibit the expression of their own genes, giving rise to oscillatory patterns in protein concentrations and hence, circadian rhythms.1,2,18 In mammals, core clock genes include Per, Cry, Bmal1 and Clock genes. Their proteins act by inhibiting or stimulating transcriptions of other core clock genes. The proteins of the Bmal1 and Clock genes, namely BMAL1 and CLOCK, form a complex CLOCKBMAL1 inside the nucleus. This complex activates the transcription of the Per and the Cry genes. In the cytoplasm, Per and Cry RNA translate to their respective proteins, PER and CRY. Some of these proteins dimer- Fig. 2: Mammalian circadian clock ize to form the complex PER-CRY which returns to the mechanism (Francis Levi, EORTC nucleus where, by binding to the CLOCK-BMAL1 com- Chronotherapy Group). plex, it prevents the further transcription of Per and Cry genes. Thus a negative feedback loop is created; PER and CRY proteins blocking the transcription of their own genes. The above mechanism for the mammalian clock is illustrated in Fig. 2. Pacific Symposium on Biocomputing 13:402-413(2008) September 26, 2007 22:2 WSPC - Proceedings Trim Size: 9in x 6in agarwal 2.2. Models of Drosophila and mammal circadian rhythms Computational models are available for circadian rhythms in Drosophila, Neurospora, cyanobacteria and mammalian systems.1,2 These models, useful for computing concentrations of core clock genes, take into account the processes of transcription, translation and phosphorylation. The mammalian circadian clock model we use1 consists of 16 variables (hence 16 differential equations) and 52 parameters. It incorporates the effect of negative autoregulation of Per/Cry gene expression by their own proteins. The Drosophila circadian model we use,2 consisting of only 5 variables (5 equations) and 18 parameters, is small enough to be reproduced here: d Mp Mp Kn = vs n I n - vm , dt KI + PN Km + Mp d P0 dt d P1 dt d P2 dt d PN dt = ks M p - v1 P0 P1 + v2 , K1 + P0 K2 + P1 P0 P1 P1 P2 = v1 - v2 - v3 + v4 , K1 + P0 K2 + P1 K3 + P1 K4 + P2 P1 P2 P2 = v3 - v4 - vd - k1 P2 + k2 PN , K3 + P1 K4 + P2 Kd + P2 = k1 P2 - k2 PN . Fig. 3: Drosophila clock mechanism2 . (1) The variables in Eq. 1 are the same as those shown in Fig. 3; parameters that lead to endogenous oscillations are taken from Gonze/Goldbetter.2 Note that Eq. 1 is in the canonical nonlinear ODE/DAE system form d q(x(t)) + f (x(t)) + b(t ) = 0, (2) dt with q(x) x and b(t ) 0. b(t ) represents the influence of external inputs, such as light, which affects the transcription rate of the Per gene in both mammals and Drosophila. The effect of light is modelled by including a parameter in the rate equations of the Per gene (vs in Eq. 1), which we recast in the form of Eq. 2 with b(t ) = 0 (details in Section 4). 3. Oscillators and PPV phase macromodels The quantitative study and design of oscillators has a rich history in engineering, particularly in electronics: oscillators are fundamental components in virtually all electronic systems. For example, they are widely used in communication systems for frequency translation of information signals; phase locked loops (PLLs) for clock generation and frequency synthesis, etc.. As noted earlier, phase macromodelling techniques are widely used to improve simulation efficiency and accuracy5,19 in electronics. In particular, the Perturbation Projection Vector (PPV) phase macromodel6,9 is well established on account of its rigorous Floquet-theoretic underpinnings, broad applicability, effective numerical extraction procedures, large simulation speedups and extensive validation. We have already noted its advantages in Section 1; here, we summarize mathematical details of the model. For expositional convenience, we assume an ODE form for an oscillator under external perturbation: d x(t) + f (x(t)) = b(t ). (3) dt b(t ) is the vector of perturbations applied to the free running oscillator; x(t) and f (x(t)) have their usual meanings, as in Eq. 2. The solution of this perturbed oscillator can be shown6 to be in the form x p (t ) = xs (t + (t )) + y(t + (t )), (4) where xs (t ) is the periodic, oscillatory solution of the unperturbed oscillator and (t ) is a phase deviation caused by the external perturbation b(t ). y(t + (t )) is an amplitude Pacific Symposium on Biocomputing 13:402-413(2008) September 26, 2007 22:2 WSPC - Proceedings Trim Size: 9in x 6in agarwal variation; it is typically very small in circadian oscillators2 and is therefore of secondary importance compared to the phase deviation (t ). Using a nonlinear extension of Floquet theory, Demir et al6 proved that (t ) is governed by the scalar, nonlinear, time-shifted differential equation (t ) = vT (t + (t )) · b(t ), (5) 1 where vT (t ) is a periodic vector known as the perturbation projection vector or PPV. Im1 portantly, they also showed that the PPV can be calculated efficiently via simple postprocessing steps following time- or frequency-domain steady-state computation.6,9 Each component of the PPV waveform represents the oscillator's "nonlinear phase sensitivity" to perturbations of that component. The PPV needs to be extracted only once from Eq. 1 (even if parameters change, see the description of parameterized PPVs below); once extracted, Eq. 5 is used for simulations. 3.1. Using the PPV macromodel for systems of coupled oscillators By employing b(t ) in Eq. 5 to capture coupling, PPV macromodels can be composed to represent systems of many coupled oscillators with different characteristics. For purposes of illustration, we outline the procedure for N identical oscillators coupling via only one component of b(t ). This results in the following set of governing equations for the coupled system: i (t ) = vT (t + i (t )) · i (t ), i 1, · · · , N , (6) where i (t ) is the phase shift of oscillator i, v(t ) is the phase sensitivity of the node on which coupling occurs and i (t ) is the perturbation resulting on oscillator i due to coupling from other oscillators. If the coupling i (t ) and phase sensitivity v(t ) are purely sinusoidal, it is easy to show that Eq. 6 is equivalent to Kuramoto's model.15 In general, however, Eq. 6 is far more accurate since it considers all harmonics of the PPVs. We use the coupling function model given in To et al20 as i (t ) in Eq. 6, and solve for the phase dynamics of a 20×20 network of coupled oscillators. 3.2. Injection locking analysis When an external signal of frequency f is injected into an oscillator with a central frequency fo close to f , the oscillator can lock to the injected signal both in phase and frequency. This phenomenon is known as injection locking and can be very easily captured by the PPV macromodel of the oscillator.4 It has been shown5 that when injection locked, an oscillator's phase shift (t ) varies linearly with time as (t ) (t ) = t+ , (7) o o where o is the natural frequency of the unperturbed oscillator and the difference between the frequencies of the injected signal and the unperturbed oscillator. (t ) represents a bounded, periodic phase difference function, the exact form of which can be determined via time-course or steady-state simulation5,8 of Eq. 5. The presence of injection locking can therefore be detected by comparing the time-average of (t ) with . o 3.3. Parameterized PPV macromodels Circadian rhythm models typically involve large numbers of model parameters. For example, there are 18 parameters in the Drosophila clock model,2 while the mammalian clock model1 has 52 parameters. The values of these parameters are chosen so that the model's predictions best fit experimental observations. Leloup/Goldbetter17 have noted that circadian rhythm properties (particularly frequency) are sensitive to variations in several parameters. The conventional approach to assessing the effect of parameter variations involves brute-force time-course simulation of circadian models, a process that is not only expensive but can also generate numerical inaccuracies in phase.5 Pacific Symposium on Biocomputing 13:402-413(2008) September 26, 2007 22:2 WSPC - Proceedings Trim Size: 9in x 6in agarwal We, instead, use an extended form of Eq. 5 that directly incorporates parameter variations ­ we call this the parameterized PPV macromodel.13 The key advantage of the parameterized PPV macromodel is that it does not involve re-extracting the PPV when parameters change - this leads to huge speedups when, e.g., many coupled oscillators with different parameters are involved. The parameterized PPV equation is given by (t ) = vT (t + (t )) · (b(t ) - S p (t + (t )) p), (8) 1 where p is a vector containing parameter variation terms and S p (t ) is a periodic, timevarying matrix function given by f S p (t ) = | (9) . p xs (t ), p repIn Eq. 9, xs (t ) denotes the natural periodic solution of the unperturbed oscillator; p resents the vector containing nominal (basal) parameter values. This extra term captures phase deviations due to parameter variations, without having to re-extract the PPV when the parameters change. It also enables the study of the effects of multiple parameters varying at the same time. 4. Simulation of mammalian and Drosophila melanogaster circadian rhythms using PPV macromodels In this section, we present results obtained by applying PPV macromodelling, described in Section 3, to mammalian and Drosophila circadian rhythm models.1,2 We first extract PPV macromodels for both circadian systems at nominal parameters9 and then simulate for phase deviation with external perturbation to demonstrate injection locking. We model the external perturbations as changes in external light intensity by first assigning a constant value to the light sensitive parameter vs (signifying darkness) and then applying an external light signal of intensity (10) L(t ) = A + Asin( t )W /m2 , where = 2 f , f being the frequency of the light/dark cycles, i.e., corresponding to 1 cycle in 24 hours. Often, light is modelled as a step function for simulations in biological systems (i.e., constant values for light and dark conditions respectively). However, to correspond more closely with continuously changing light intensities in reality and to illustrate the generality of the PPV model, we apply sinusoidal intensity waveform around an average value.21 (Note that any other shape, including step function shapes, can be handled equally easily). We assume the experimental setup used by Usui/Okazaki,21 where the illuminance of light is varied from 20 lux to 0.01 lux (i.e., variation in light intensity from 0.15 W /m2 to 0.00009 W /m2 ), giving A 0.05W /m2 in Eq. 10. Moreover, Eq. 10 multiplied by a constant gives the term b(t ) of Eq. 2, where the constant signifies the change in Per gene concentration for 1 W /m2 of light intensity. In this paper, we assume the constant to be equal to 1 nM /(W /m2 ). The constant can be modelled accurately in future experiments. We also extract parameterized PPV macromodels to study the effect of parameter variations in two cases ­ with and without external light variations. In the absence of external light variations, phase deviations from the parameterized PPV macromodel are useful for predicting changes in free-running frequency. When external light perturbations are included in parameter-varying PPV simulations, lock range information is also generated. Finally, we put the above single-oscillator PPV macromodels together to model a locally-coupled 20×20 network of oscillators ­ a simple representation of a spatially multiscale, coupled circadian system. We use this model to demonstrate synchronization behavior, obtaining speedups of about 240× over traditional time-course simulation. 4.1. Time-course simulations using full ODE models For reference and validation, we first perform time-course simulations of the two ODE circadian rhythm models directly, to obtain concentration waveforms for all clock proteins and mRNAs in the model. The waveforms thus obtained are shown in Fig. 4(a) and Fig. 4(b). We observe an anti-phase relationship between the concentrations of the Per/Cry and Bmal1 mRNAs, as expected from theory.1 The period of the oscillating waveforms is equal to 23.8 hrs for the mammalian clock and 22.4 hrs for the Drosophila clock. Pacific Symposium on Biocomputing 13:402-413(2008) September 26, 2007 22:2 WSPC - Proceedings Trim Size: 9in x 6in agarwal Circadian Oscillations (mammals) Concentration (nM) 8 6 4 2 0 450 500 550 time (hrs) 600 Concentration (nM) Circadian Oscillations (Drosophila) 1.5 1 0.5 Per Cry Bmal1 260 280 300 320 340 360 380 time (hrs) (a) Per, Cry and Bmal1 Gene Concentra- (b) Per tions (Mammals) (Drosophila) Gene Concentration Fig. 4. (a) Plot of core clock gene concentrations (Per, Cry and Bmal1) in mammals vs. time. The concentrations are oscillatory and there is an antiphase relationship between the Per/Cry and Bmal1 concentrations. (b) Plot of the Per gene concentration in Drosophila vs. time 4.2. Circadian PPV macromodels In this section, we extract the PPV macromodel of the circadian oscillator for both models. Fig. 5(a) and Fig. 5(b) show the PPV waveforms of Per gene concentrations. This waveform gives the phase sensitivity of the concentration at each time instant and can be directly used to find the new concentration waveform under the effect of an external perturbation. It is equivalent to the phase response curve described by Winfree,14 with the only exception that PPV waveforms do not involve sinusoidal simplifications,15 implying greater accuracy, as already noted previously. By inspecting the phase sensitivity at each time instant, it becomes possible to determine the time at which light should be applied to shift the oscillator's time-keeping forward or backward. At zero crossings of the PPV phase sensitivity function, for example, a light pulse will have no effect on the phase/frequency characteristics of the oscillator. PPV waveform (Per gene, mammals) Phase sensitivity Phase sensitivity 5 10 15 time (hrs) 20 PPV waveform (Per gene, Drosophila) 10 5 0 -5 -10 0 5 10 15 time (hrs) 20 2 1 0 -1 0 (a) PPV waveform for Per mRNA (gene) (b) PPV waveform for Per mRNA (gene) concentration (Mammals) Concentration (Drosophila) Fig. 5. (a) and (b) Plot of PPV phase sensitivities vs. time for the Per gene concentration. Speedups: (a) For the mammalian clock model: full time course simulation using the Backward Euler (BE) integration method requires 18 seconds of computer time, while finding the free-running steady state via harmonic balance analysis takes about 6 seconds. This is followed by the PPV extraction algorithm, which takes around 1.5 seconds; the total time required for PPV extraction is about 7.5 seconds, representing a speedup of some 2.5 ×. (b) For the Drosophila model: full time course simulation takes about 13 secs; harmonic balance analysis takes 4 seconds, PPV extraction 0.5 seconds; resulting in a speedup of 3×. To gauge the accuracy of the PPV macromodels, we plot concentration waveforms for the Per gene, obtained from time course and PPV simulations, in Fig. 6. Fig. 6(a) shows waveforms for an locked oscillator (distinct frequencies), while Fig. 6(b) shows waveforms Pacific Symposium on Biocomputing 13:402-413(2008) September 26, 2007 22:2 WSPC - Proceedings Trim Size: 9in x 6in agarwal for a locked oscillator. (a) (b) Fig. 6. Plots of Per gene (mRNA) concentration obtained from transient and PPV macromodel simulations. (a) Unlocked case. (b) Locked case. 4.3. Simulation of injection locking In order to study the effects of external perturbations on circadian rhythms, we calculate phase deviations due to external perturbations (Eq. 10) by solving Eq. 7. If the period of the externally applied signal is close enough to the oscillator's free-running frequency, entrainment or injection locking occurs. 4.3.1. Mammalian clock model The free-running frequency of the mammalian circadian clock is f o = 4.19 × 10-2 hr-1 . We apply an external signal with of = -0.00623. Fig. 7(a) depicts injection locking with f a light input of 0.009 + 0.009sin( t ) W /m2 , where the locking starts around 690 hours ( 30 cycles). Fig. 7(b) shows the same curve for light input of 0.05 + 0.05sin( t ) W /m2 ; locking time reduces to about 260 hrs ( 11 cycles). From these results, we can infer that with smaller light intensities, the resetting phenomenon takes a longer time. Phase deviation (slope=-0.006993) Phase deviation (hrs) 0 -20 -40 -60 2000 4000 6000 8000 10000 time (hrs) Phase deviation (hrs) Phase deviation (slope=-0.007000) 0 -20 -40 -60 2000 4000 6000 8000 10000 time (hrs) (a) t vs time (b) t vs time Fig. 7. (a) Plot of phase deviation (t ) vs. time for the mammalian clock model, with light input 0.009 + 0.009 sin( t )W /m2 . The slope of -0.0069 indicates injection locking, with lock reached in about 690 hours. (b) Plot of (t ) vs. time for the mammalian clock model, with light input 0.05 + 0.05 sin( t )W /m2 . The slope is -0.007; lock is reached in about 260 hours. 4.3.2. Drosophila clock model The free-running frequency of the Drosophila circadian oscillator was f o = 4.48x10-2 hr-1 ; the frequency of the injected light signal was f = 4.16x10-2 hr-1 , leading to of = -0.071. The intensity of the applied light is given by Eq. 10, with A = 0.05W /m2 . f Fig. 8(a) shows the phase deviation vs. time; its slope is 0.071, equal to the relative frequency difference, thus confirming that the oscillator is locked to the injection frequency. Pacific Symposium on Biocomputing 13:402-413(2008) September 26, 2007 22:2 WSPC - Proceedings Trim Size: 9in x 6in agarwal Frequency Difference (%) Phase deviation (hrs) Phase deviation (slope=-0.071000) 0 -200 -400 -600 2000 4000 6000 8000 10000 time (hrs) Freq. Difference between the oscillator external input 0.2 0.1 0 2000 4000 6000 8000 10000 time (hrs) (a) t vs time (b) Frequency Deviation vs time Fig. 8. (a) Plot of phase deviation (t ) vs. time for the Drosophila clock model. The slope is -0.071 (injection locked) (b) Plot of frequency deviation vs. time. 4.4. Lock range vs. injection amplitude We also calculate the lock range (frequency range Locking Range vs Injection Amplitude over which the oscillator remains locked to the external 0.6 signal) for the mammalian clock and plot it as a function of injection amplitude. We find that the lock range 0.4 increases roughly linearly with injection amplitude, as can be seen in Fig. 9. However, at higher amplitudes, 0.2 the linearity between lock range and injection amplitude collapses. By calculating the lock range for a given 0.2 0.4 0.6 light amplitude, we can infer whether the system would Injection Amplitude (W/m2) lose its rhythmicity or not on exposure to that particular g 9 L o king Ra vI light. Conversely, one can calculate the light amplitude Fim.pl:tudec(Mammngie n sClnjection A i al a ock required to synchronize the free running oscillators. Speedups: (a) Mammalian clock model: Time course simulations take 18 seconds. PPV macromodel simulations take 2 seconds after PPV extraction, resulting in a speedup of 9×. (b) Drosophila clock model: Time course simulations take 13 seconds; PPV macromodel simulations take 1 second; resulting in a speedup of 13×. Relative frequency difference 4.5. Parameter variation simulations To study the effect of parameter variations on circadian rhythms, we first simulate the phase deviations given by Eq. 8 with b(t ) = 0. As an example, we vary all parameter values by 10% of their nominal values; i.e., p = 0.1 p. The slope of the phase deviation curve gives the relative change in frequency due to the change in parameters. As is evident from Eq. 8, we can study the effects of all possible combinations of parameter variations. For the mammalian clock model, the relative frequency change is found to equal 0.186; i.e., the new frequency is equal to 20.1 hours (Fig. 10(a)). For the Drosophila clock model, the relative frequency change equals -0.114, i.e., the new frequency equals 25.2 hrs (Fig. 10(a)). It is evident that even small changes in parameter values can affect rhythm frequency significantly. Next, we combine parameter variations with external perturbations to the oscillator. Using the slope of the phase deviation curve, we can calculate the range over which parameters can be varied while keeping the oscillator oscillator still locked to the injection signal frequency for the same light input. As an example, we vary all parameters simultaneously and find the respective variation range for each model. In case of the Drosophila model, parameters can be varied from -5t o10% without loss of lock; while for the mammalian model, the range is smaller, -5t o1% variation. The light input in both the cases is given by Eq. 10 with A = 0.05W /m2 . Pacific Symposium on Biocomputing 13:402-413(2008) September 26, 2007 22:2 WSPC - Proceedings Trim Size: 9in x 6in agarwal Phase deviation (slope=0.186167) Phase deviation (hrs) 500 400 300 200 100 0 500 1000 1500 2000 2500 time (hrs) Phase deviation (hrs) Phase deviation (slope=-0.114239) 0 -100 -200 -300 500 1000 1500 2000 2500 time (hrs) (a) t vs time (b) t vs time Fig. 10. (a) Plot of the phase deviation (t ) vs. time for the mammalian clock model with 10% variation in the parameter values. The slope of the curve = 0.186 implying that the new frequency of oscillations equals to 20.1 hrs. (b) Plot between the phase deviation (t ) and time for the Drosophila Clock Model with 10% variation in the parameter values. The slope = -0.114 and the new frequency of oscillations is hence equal to 25.2 hrs. 4.6. Synchronization of coupled oscillators In this section, we extend the single oscillator analysis to a system of many coupled oscillators (i.e., a system of several interacting biological cells, each behaving as an individual oscillator and oscillating with a period of 24 hrs). We consider a system of 400 mammalian clock oscillators arranged in a 20x20 grid, as shown in Fig. 11.1 The oscillators are identical in all respects except for their free-running frequencies, which are selected randomly from a uniform distribution. Each oscillator is modelled by a system of 16 ODEs (as used before for single oscillator analyses). In order to introduce coupling between the oscillators, we use a recently proposed coupling model given in To et al,20 wherein neurotransmitters act as synchro- Fig. 11: 2-dimensional oscillator grid. The nizing agents between the cells. Then, using a numbers indicate the weight factors used for PPV macromodel for each oscillator augmented the coupling. The black solid circle represents by coupling equations, we simulate the entire a particular cell of interest.20 oscillator system. We use Eq. 6 to calculate the phase deviations for each oscillator, recording instantaneous phases at regular intervals. In every phase plot (e.g., as shown in Fig. 12(a)), a small rectangle represents an individual oscillator; the colour of the rectangle represents its phase visually; e.g., dark red denotes a phase of , while dark blue denotes 0 phase. Fig. 12(a) and Fig. 12(b) show phase plots at t = 1T and 5T respectively (T is the free running frequency of an oscillator) in the absence of coupling. The absence of coupling can easily be surmised, from the random nature of the plots (absence of any pattern, i.e., unsynchronized phases). For the coupled case, Fig. 12(d) shows the phases at 0.5T , when all the oscillators start synchronizing to the same phase (and frequency). Fig. 12(e) and Fig. 12(f) are the phase plots at later stages, confirming synchronization amongst the coupled oscillators. We have also varied the random center frequency distributions of the oscillators, and found that with the same coupling strength, the oscillators cease to lock to each other for deviations greater than 0.5 of the free-running circadian period Speedups: (a) Time course simulations require 12 hours for full simulations, including the time required for the formation of the coupling matrix. (b) PPV simulations require 158 seconds for complete simulations. Hence, we obtain a speedup of 240×. If the system size is larger and the oscillator model is more complex, the speedups will be greater. Pacific Symposium on Biocomputing 13:402-413(2008) September 26, 2007 22:2 WSPC - Proceedings Trim Size: 9in x 6in agarwal 2 2 2 4 4 4 6 6 6 8 8 8 10 10 10 12 12 12 14 14 14 16 16 16 18 18 18 20 2 4 6 8 10 12 14 16 18 20 20 2 4 6 8 10 12 14 16 18 20 20 2 4 6 8 10 12 14 16 18 20 (a) Phase Plot at t = 1T (No Coupling) 2 2 (b) Phase Plot at t = 5T (No Coupling) 2 (c) Phase Plot at t = 0 (Coupled) 4 4 4 6 6 6 8 8 8 10 10 10 12 12 12 14 14 14 16 16 16 18 18 18 20 2 4 6 8 10 12 14 16 18 20 20 2 4 6 8 10 12 14 16 18 20 20 2 4 6 8 10 12 14 16 18 20 (d) Phase Plot at t = 0.5T (e) Phase Plot at t = 0.75T (f) Phase Plot at t = 1.25T (Coupled) (Coupled) (Coupled) Fig. 12. (a) and (b) Phase plots in case of no intercellular coupling between individual oscillators. (c) - (f) Phase plots showing the synchronization of coupled oscillators (all oscillators at the same phase) 5. Conclusion We have applied PPV phase macromodelling techniques to mammalian and Drosophila circadian rhythms, for the first time. These techniques provide fast/accurate simulations of oscillator systems, predicting synchronization and resetting in circadian rhythms via injection locking cued by light inputs. In addition, PPV waveforms provide direct insight into the effect of light on phases of the oscillating rhythms. We have accurately predicted synchronization in a coupled multi-scale system of 400 circadian oscillators using PPV macromodels. Finally, the efficacy of parameterized PPV macromodels for circadian problems has also been demonstrated. References 1. J.C. Leloup and A. Goldbetter, "Toward a detailed computational model for the mammalian circadian clock" in Proceedings of the National Academy of Sciences of the United States of America , 7051(June 2003). 2. J. D. Gonze and A. Goldbetter, "Robustness of circadian rhythms with respect to molecular noise" in Proceedings of the National Academy of Sciences of the United States of America , 673(January 2002). 3. Y. Touitou, Biological Clocks: Mechanisms and Applications (Proceedings of the International Congress on Chronobiology) (Elsevier, Paris, France, 1997). 4. R. Adler, "A study of locking phenomenon in oscillators" in Proceedings of the I.R.E. and Waves and Electrons 34, 351 (1946). 5. X. Lai and J. Roychowdhury, "Capturing Oscillator Injection Locking via Nonlinear Phase-Domain Macromodels" in IEEE Trans. MTT 52, 2251(September 2004). 6. A. Demir, A. Mehrotra and J. Roychowdhury, "Phase noise in oscillators: a unifying theory and numerical methods for characterization" in IEEE Trans. Ckts. Syst. ­ I: Fund. Th. Appl. 47, 655(May 2000). 7. T. Mei and J. Roychowdhury, "A Robust Envelope Following Method Applicalbe to both Non-autonomous and Oscillatory Circuits" in(July 2006). 8. T. Mei and J. Roychowdhury, "An Efficient and Robust Technique for Tracking Amplitude and Frequency Envelopes in Oscillators" in(November 2005). 9. A. Demir and J. Roychowdhury, "A Reliable and Efficient Procedure for Oscilla- Pacific Symposium on Biocomputing 13:402-413(2008) September 26, 2007 22:2 WSPC - Proceedings Trim Size: 9in x 6in agarwal 10. 11. 12. 13. 14. 15. 16. 17. 18. 19. 20. 21. tor PPV Computation, with Phase Noise Macromodelling Applications" in IEEE Trans. Ckts. Syst. ­ I: Fund. Th. Appl. , 188(February 2003). X. Lai and J. Roychowdhury, "Fast, accurate prediction of PLL jitter induced by power grid noise" in(May 2004). X. Lai and J. Roychowdhury, "Fast Simulations to Large Networks of Nanotechnological and Biochemical Oscillators for Investigating Self-Organization Phenomenon" in Proc. IEEE ASP-DAC (2006). A. Demir and J. Roychowdhury, "A reliable and efficient procedure for oscillator PPV computation, with phase noise macromodelling applications." in IEEE Transactions on Computer-Aided Design 22, 188(February 2003). X. Z. Wang and J. Roychowdhury, "PV-PPV: Parameter Variability Aware, Automatically Extracted, Nonlinear Time-Shifted Oscillator Macromodels" in(June 2007). A. Winfree, "Biological Rhythms and the Behavior of Populations of Coupled Oscillators" in Theoretical Biology 16, 15 (1967). Y. Kuramoto, "Chemical osillations, waves and turbulence" in Springer (1984). S.H. Strogatz, "From Kuramoto to Crawford: exploring the onset of synchronization in populations of coupled oscillators" in Elsevier 143, 1 (2000). J.C. Leloup and A. Goldbetter, "Modeling the mammalian circadian clock: Sensitivity analysis and multiplicity of oscillatory mechanisms" in Theoretical Biology , 5 4 1 ( Ap r i l 2 0 0 4 ) . Circadian rhythms laboratory homepage http://www.circadian.org/main.html. X. Lai and J. Roychowdhury, "Macromodelling Oscillators Using Krylov-Subspace Methods" in Proc. IEEE ASP-DAC (January 2006). T.L. To, M.A. Henson, E.D. Herzog and F.J. Doyle III, "A Molecular Model for Intercellular Synchronization in the Mammalian Circadian Clock" in Biophysical Journal 92, 3792 (2007). Y. S. Usui and T. Okazaki, "Range of entrainment of rat circadian rhythms to sinusoidal light-intensity cycles" in AJP-Regulatory, Integrative and Comparitive Physiology 278, R1148(May 2000).