Dianailys Nu Ñez-Reyes

Jacek Kłos

Millard H Alexander

Paul J Dagdigian
email: pjdagdigian@jhu.edu

Kevin M Hickson
email: kevin.hickson@u-bordeaux.fr

Experimental and theoretical investigation of the temperature dependent electronic quenching of O( 1 D) atoms in collisions with Kr

The kinetics and dynamics of the collisional electronic quenching of O( 1 D) atoms by Kr has been investigated in a joint experimental and theoretical study. The kinetics of quenching was measured over the temperature range 50 -296 K using the Laval nozzle method.

O( 1 D) atoms were prepared by 266 nm photolysis of ozone, and the decay of the O( 1 D) concentration was monitored through vacuum ultraviolet fluorescence at 115.215 nm, from which the rate constant was determined. To interpret the experiments, a quantum closecoupling treatment of the quenching transition from the 1 D state to the 3 P j fine-structure levels in collisions with Kr, and also Ar and Xe, was carried out. The relevant potential energy curves and spin-orbit coupling matrix elements were obtained in electronic structure calculations. We find reasonable agreement between computed temperature-dependent O( 1 D)-Rg (Rg = Ar, Kr, Xe) quenching rate constants and the present measurements for Kr and earlier measurements. In particular, the temperature dependence is well described.

I. INTRODUCTION

The collisional de-activation of electronically excited O( 1 D) atoms is an important process in the chemistry of the earth's atmosphere. The ultraviolet photolysis of ozone in the stratosphere and troposphere leads to the production of oxygen atoms in this excited state, as well as in the ground 3 P state. The relatively low steady-state concentration of O( 1 D) atoms in the atmosphere is controlled by the rate of photolytic production and loss by collisional electronic quenching.

While not as relevant for atmospheric chemistry as de-activation by and chemical reaction with molecular collision partners, the quenching of O( 1 D) by rare gas (Rg) atoms is a relatively simple collisional process and has attracted much interest:

O( 1 D) + Rg → O( 3 P j ) + Rg (1)
In Eq. ( 1), the subscript j denotes that the quenched O( 3 P) atoms can be formed in the finestructure levels with angular momenta j = 2, 1, and 0.

The rate constant for the electronic quenching of O( 1 D) by the light rare gases is known to be small. 1 There have been a number of measurements of the O( 1 D) quenching rate constant by the heavier rare gases Ar (Refs. Similarly, there have been measurements of the electronic quenching of the isovalent S( 1 D)

atom by Ar atoms 8 and other species. 9,10 In the recent work by Lara et al., 8 the rate constant for quenching of S( 1 D) by Ar was measured as a function of temperature from 298 down to 5.8 K. This rate constant was compared with a theoretical rate constant computed through quantum scattering calculations using calculated S( 1 D, 3 P) potential energy curves and spin-orbit matrix elements.

Previous measurements of the rate constant for electronic quenching of O( 1 D) by Kr have been reported down to a temperature of 113 K. 2 In the present study, we report measurement of this rate constant from 50 to 296 K. As a compliment to our measurement and previous determinations of the rate constants for the quenching of O( 1 D) by the heavy rare gases, we have also computed integral cross sections and rate constants for quenching by Ar, Kr, and Xe. We employed the previously computed 7 O( 1 D, 3 P) potential energy curves and spin-orbit matrix elements and have carried out here similar calculations for O( 1 D, 3 P)-Kr and O( 1 D, 3 P)-Ar. We find that the energy-dependent cross sections for all three collision pairs show many sharp peaks, that are most prominent for Kr and Ar as the collision partner. We show that these peaks are resonances that are quasi-bound by the centrifugal barrier at the highest partial wave that contributes to the quenching cross section.

This paper is organized in the following manner. Section II presents the experimental method

II. EXPERIMENT

A. Methods

The experimental part of this work was performed using a supersonic flow (Laval nozzle) reactor operating in continuous mode. A detailed description of the apparatus can be found in earlier work. 12,13 Three different axisymmetric nozzles were employed during this investigation, allowing generation of flows at three different fixed low temperatures (50 K, 75 K, and 127 K), each with a uniform density profile and a constant velocity. Although we have a range of Laval nozzles available with Ar, N 2 and N 2 /SF 6 mixtures as the carrier gases, only Ar based nozzles could be used in this study due to the efficient non-reactive quenching of O( 1 D) by N 2 below room temperature. 5,14 The supersonic flow characteristics of these nozzles are presented in Table 1 of Grondin et al. 5 and were obtained in separate calibration experiments. In these experiments, the Mach number (the ratio of the flow velocity to the speed of sound for the specific carrier gas at the local temperature) was determined from standard isentropic flow equations relating the stagnation pressure (reservoir pressure) and the impact pressure obtained with a Pitot tube inserted into the supersonic flow. Measurements were also performed at room temperature by operating without a nozzle and by throttling the pumps to reduce the flow velocity, thereby eliminating temperature and pressure gradients within the reactor. to generate tuneable VUV light, as well as to detect resonant VUV emission, are described in previous work. 5,16,17 The O( 1 D) VUV LIF signal was followed as a function of delay time between photolysis and probe lasers, with a boxcar integration system coupled to a PC used for signal acquisition and processing. Both lasers and the acquisition electronics were synchronized using a digital delay generator. A total of 30 laser shots were recorded at each time interval, with a least 70 time intervals used to establish the O( 1 D) decay profiles. In addition, at least 15 negative time points (that is with the probe laser firing prior to the photolysis laser) were recorded for each decay profile to set the baseline level. The gases O 2 (99.999%), Ar (99.999%), Kr (99.999%) and Xe (99.998%) were all used with no further purification during the present experiments. The co-reagent Kr concentrations were determined from the product of its flow ratio (F Kr / F tot ) and the supersonic flow density. Calibrated mass flow controllers were used to drive the gas flows, with calibration factors determined by a pressure rise at constant volume method for the specific gas used.

B. Results

The concentration of Kr used in the experiments was always much greater than the O( 1 D) concentration ([O 2 ] max = 7.0×10 13 cm -3 with [Kr] min = 1.3×10 15 cm -3 ) so that the pseudo-firstorder approximation was valid at all temperatures. Under these conditions the O( 1 D) VUV LIF signal decreased exponentially as a function of time according to the formula I t = I 0 exp(-k ′ t), where I t and I 0 are the time dependent and initial O( 1 D) VUV LIF signals, respectively, and t is the time. The quantity k ′ is the pseudo-first-order loss rate for O( 1 D) atoms and contains contributions from several different processes: 5 When Kr is added to the system ([Kr] = 1.48×10 16 cm -3 for the blue data points in Fig. 1), O( 1 D) atoms decay much more rapidly, with a pseudo first-order rate constant

O( 1 D) + Kr → O( 3 P) + Kr (2) O( 1 D) + Ar → O( 3 P) + Ar (3) O( 1 D) + O 2 → O( 3 P) + O 2 (4) O( 1 D) + O 3 → products (5) O( 1 D) → diffusion (6) such that k ′ = k O( 1 D)-Kr [Kr] + k O( 1 D)-Ar [Ar] + k O( 1 D)-O 2 [O 2 ] + k O( 1 D)-O 3 [O 3 ] + k diff .
k ′ = (1.78 ± 0.8) × 10 5 s -1 .
Pseudo-first-order rate constants were determined over a wide range of different Kr concentrations at each temperature. A plot of the pseudo-first-order rate constants against the corresponding Kr concentrations yielded the second-order rate constant from a weighted linear least-squares fit to the data. The second order plots obtained at 50 K and at 296 K are shown in Fig. 2.

In these experiments, the fits were weighted by the uncertainty derived by exponential fits to first-order decay profiles such as those shown in Fig. 1. The large y-axis intercepts of these secondorder plots represent O( 1 D) loss essentially through reaction (3), as described above. In addition, since large Kr flows were required in some of these experiments due to the relatively slow O( 1 D) quenching rate, it was necessary to reduce the Ar flow rate to maintain a constant reservoir pressure (a prerequisite for experiments of this type for the Laval nozzle calibrations to be considered valid).

Consequently, as [Ar] decreased slightly with increasing [Kr], a small correction was applied to the O( 1 D) first-order decay rates to compensate for the reduced relaxation by Ar. This correction was given by the following formula:

([Ar] 0 -[Ar] Kr ) k O( 1 D)-Ar
, where [Ar] 0 is the nominal Ar density for the specific nozzle, [Ar] Kr is the Ar density with added Kr, and k O( 1 D)-Ar is the measured secondorder rate constant for O( 1 D) removal by Ar. 2,5 Adding this correction to the measured pseudofirst-order rate constants resulted in an increase of less than 10% to the second-order rate constants.

The corrected second-order rate constants obtained for the O( 1 D) + Kr reaction are presented in Table I, with other relevant information, and are displayed as a function of temperature in Fig. 3 along with earlier results.

III. THEORY A. Electronic structure calculations

In order to perform quenching scattering calculations, we need information on the potential energy curves (PEC's) and spin-orbit matrix elements for the interaction of closed-shell 1 S rare gas atoms with the oxygen atom in its ground 3 P and first excited 1 D electronic states. We employ similar computational methodology to that described in earlier work on O( 1 D, 3 P)-Xe quench- interaction is described by five PEC's (three unique) belonging to 1 Σ + , 1 Π (doubly degenerate) and 1 ∆ (doubly degenerate) representations. We use the MOLPRO suite of programs [START_REF] Werner | Molpro, version 2012.1, a package of ab initio programs[END_REF] to obtain all of these PEC's from first principles calculations.

As in our work on O-Xe quenching, 7 we perform two sets of ab initio calculations of the O-Kr and O-Ar PEC's. The first set is based on the multi-reference configuration interaction with single and double excitations [START_REF] Werner | [END_REF]20 method with a correction for quadruple excitations 21

(MRCISD+Q).

The MRCISD+Q calculations also provide necessary spin-orbit coupling terms between all the electronic states considered in the O-Rg systems. The second set employs the spin unrestricted coupled-cluster method with single, double and non-iterative triple excitations [UCCSD(T)]. The UCCSD(T) method describes well all the potential energy curves except the one for the 1 Π state. 7 For that state we will use the PEC obtained from the set of MRCISD+Q calculations.

We perform all-electron ab initio calculations. In these, we have to account for the scalar relativistic effects that can be important for the Kr atom. For this reason we employ the aug-mented correlation-consistent polarized quadruple and quintuple zeta atomic basis sets designed for Douglas-Kroll relativistic calculations [aug-cc-pvqz-DK in the MRCISD and aug-cc-pv5z-DK in the UCCSD(T) calculations].

The MRCISD calculations were started from the state-averaged multiconfiguration self consistent field calculations (MCSCF). The active space employed in the MCSCF calculations included 15 core orbitals, 5 closed (doubly occupied), but correlated, orbitals, and three active p orbitals of the oxygen atom. The states of Σ, Π and ∆ symmetry were obtained by applying the projection operator methodology implemented in the MOLPRO program, ensuring that the wavefunctions are eigenvectors of the L 2 operator.

To determine the UCCSD(T) interaction energies that describe the PEC's, we use the supermolecular method in which we correct for the basis set superposition error by applying the counterpoise procedure of Boys and Bernardi 22 (all total energies calculated in the dimer basis set).

In the case of the MRCI calculations, the lack of size-consistency of the method is remedied by subtracting the interaction energy at a large (R = 30a 0 ) O-Rg separation . The correction for quadruple excitations in the MRCI calculations accounts, partially, for higher excitations, and helps to reduce the size-consistency error. II. We compare the values of D e and R e the triplet states calculated in this work with empirically determined values of these parameters derived by Aquilanti et al. 23 The latter were obtained in molecular beam experiments employing magnetically state-selected O( 3 P) atoms. The agreement of the computed and experimental values is quite good.

There have been a number of previous calculations on the triplet states of O-Kr and O-Ar. [24][25][26][27][28][29][30] The most recent of these calculations on O-Kr, which used the RCCSD(T) method with a variety of basis sets obtained the values R e = 7.417a 0 and D e = 58.3 cm -1 for the 3 Σ -state and R e = 4. Potential energy curves for the interaction of O( 3 P, 1 D) with Kr (left-hand panel) and Ar (righthand panel) in the absence of spin-orbit coupling: UCCSD(T) calculations, except for the MRCISD+Q calculation for the 1 Π state. The curves of states of singlet and triplet spin multiplicity are plotted in red and blue, respectively. The predominant spin-orbit coupling occurs between the 1 Σ + and 3 Π Ω=0 + curves, designated by the small magenta circle. An additional crossing between the 1 Σ + and 3 Σ -curves is marked with a "B". There is one additional crossing [ 1 Π -3 Σ -], but this occurs above the plotted energies. 6.516a 0 and D e = 100.7 cm -1 for the 3 Π state. Our values are in good agreement with these results.

∆ 1 Π 1 Σ + 3 Π 3 Σ - B 1 ∆ 1 Π 1 Σ + 3 Π 3 Σ - B FIG.
As noted previously, 7 there are nine and five spin projection states for the oxygen 3 P and 1 D states, respectively. The matrix elements of the 14 × 14 spin-orbit matrix expressed in a Cartesian basis were calculated at the MRCISD level using the MOLPRO suite of programs. [START_REF] Werner | Molpro, version 2012.1, a package of ab initio programs[END_REF]31 There are two unique matrix elements describing the spin-orbit coupling within the 3 P state and five such matrix elements describing the spin-orbit coupling between the 1 D and 3 P states. These matrix elements, which are defined in Ref. 7, for O-Kr and O-Ar are plotted in Fig. 5. We see in Fig. 5 that the matrix elements deviate significantly from their asymptotic values when the internuclear separation R becomes small. A similar effect was seen for O-Xe; however, the deviation from the asymptotic values was much larger. We do see in Fig. 5 that the spin-orbit matrix elements for small R are larger for Kr than for Ar. We certainly expect that the spin-orbit matrix elements would follow in the order O-Ar < O-Kr < O-Xe, going down the column of the rare gases in the periodic table.

The quenching of O( 1 D) by a rare gas can, in principle, be mediated by three different singlettriplet crossings, namely between (1) the 1 Σ + and 3 Π states, (2) the the 1 Σ + and 3 Σ -states, and

the 1 Π and 3 Σ -states. Based on our previous work on O-Xe, 7 we expect that the first crossing, which lies considerably below the O( 1 D) + Rg asymptote (denoted with small magenta circles in Fig. 4) will dominate the quenching by Kr and Ar. The second crossing also lies below the singlet asymptote and can also contribute to the quenching. In O-Xe, the third crossing lies ∼ 2000 cm -1

above the singlet asymptote and contributes, in this system, to quenching only at high collision energies. 7 . By contrast for O-Kr and O-Ar, this crossing lies much higher above the singlet asymptote, out of the range of the energies considered here (plotted in Fig. 4).

B. Scattering formalism

The formalism for quantum scattering to describe the collision-induced 1 D → 3 P transition in the oxygen atom has been presented previously. 7,32 We provide a brief description here. The Hamiltonian for the O-Rg system can be written in atomic units ( = 1) as Matrix elements between 3 P and 1 D projection states:

Ĥ(R, q) = - 1 2µR 2 d 2 dR 2 + l 2 2µR 2 + Ĥel (q) + V(R, q) + ĤSO (R, q) (7)
B yx ≡ 3 Π y |L z S z | 1 Π x , B ss ≡ 3 Σ|L z S z | 1 Σ , B xs ≡ 3 Π x |L y S y | 3 Σ , B sy ≡ 3 Σ|L y S y | 1 Π y , and B xd ≡ 3 Π x |L x S x | 1 ∆ .
where R is the vector separating the atoms, q represents the electronic coordinates, µ is the collision reduced mass, and l is the orbital angular momentum. The term Ĥel (q) is the energy of the isolated O( 3 P, 1 D) atom, and V(R, q) is the electrostatic interaction between the atoms; the latter approaches zero as R becomes larger. The term ĤSO (R, q) describes the spin-orbit interaction in the O atom; this term is affected by the approach of the rare gas atom.

The electronic part of the scattering wave function at large R is described by the quantum numbers L, S , and j, namely the electronic, spin, and total angular momentum, respectively, of the oxygen atom. We express the scattering wave function in the space frame (SF) as

|LS jlJM = m j m l jm j lm l | jlJM |LS jm j |lm l (8)
where J and M are the total angular momentum and its space-fixed projection, respectively, and l is the orbital angular momentum.

The matrix elements of V(R, q) and ĤSO (R, q) are computed in the body frame, as described in Sec. III A. The body-frame (BF) basis functions can be transformed to the SF basis with the following transformation, which has been described in detail previously: 33

|LS jlJM = [l][J] 4π 1/2 (-1) J Ω (-1) Ω          l j J 0 Ω -Ω          D J * MΩ (φ, θ, 0)|LS jΩ (9)
Here, [x] = 2x + 1, the D is a rotation matrix element, [START_REF] Zare | Angular Momentum[END_REF] , and Ω is the body-frame projection of J and j. The spin-orbit operator mixes the Ω = ±2 components of the isolated 3 P and 1 D states; this mixing must be taken into account in describing the asymptotic atomic states. We compute the spin-orbit interaction in Eq. ( 7) as ĤSO (R, q) -ĤSO (R = ∞, q) so that this interaction goes to zero for large R.

The matrix elements of V el (R, q) ≡ V(R, q) + ĤSO (R, q) are diagonal in Ω in the BF basis. Using Eq. ( 9), we can express the SF matrix elements of V el (R, q) in the following way:

L ′ S ′ j ′ l ′ JM|V el (R, q)|LS jlJM = ([l ′ ][l]) 1/2 × Ω          l ′ j ′ J 0 Ω -Ω                   l j J 0 Ω -Ω          L ′ S ′ j ′ Ω|V el (R, q)|LS jΩ (10)
The BF matrix elements of V(R, q) for the 3 P state are given in Ref. 33; these matrix elements for the 1 D state equal

LS Ω|V(R, q)|LS Ω = V Λ=Ω (R) (11)
The BF matrix elements of ĤSO (R, q) are given in the Appendix of Ref. 7.

Time-independent quantum scattering calculations for the O-Rg systems were performed with the HIBRIDON suite of programs. [START_REF] Alexander | HIBRIDON is a package of programs for the time-independent quantum treatment of inelastic collisions and photodissociation written[END_REF] Convergence of the cross sections was checked with respect to the sector width in the radial integration and the number of partial waves. Partial waves up to J = 280 were included in the calculations. The cross sections were computed on a 1 cm -1 spacing for total energies less than 2000 cm -1 and a spacing of 2 cm -1 for energies between 2000 and 5000 cm -1 .

C. Calculations The sharp features in the energy dependent cross sections plotted in Fig. 6 are intriguing. Similar sharp features were observed in computed energy dependent cross section for intramultiplet transitions in collisions of O( 3 P j ) fine-structure levels with rare gases. 29 Figure 7 displays a plot, over a narrower energy range, of the energy dependent cross sections, summed over the O( 3 P j ) fine-structure levels, for the quenching of O( 1 D) by Kr. Since the sharp features appear for each final level, for clarity we plot only the cross section summed over the fine-structure levels. We see that there are several lines of features, with the widths within a given series increasing as a function of the collision energy.

Cross Sections

To gain insight into the origin of these features, it is convenient to consider the simpler 2-state model that we employed to understand the origin of oscillations in the differential cross section for the electronic quenching of O( 1 D) by Xe. 7 As in the case of O( 3 P)-Xe, the dominant coupling enabling the quenching is between the 1 Σ + state and the Ω = 0 + component of the 3 Π state. We consider just these two diabatic potential energy curves, along with a Gaussian coupling (in units of cm -1 ):

V 12 (R) = 103 exp [-α(R -R c ) 2 ] ( 12
)
where R c = 3.825a 0 is the crossing between the two diabatic curves (circled in magenta in the left-hand panel of Fig. 4), V 12 (R) at R = R c was set to the value of the spin-orbit coupling between the two curves, and the width parameter α was set to 1.0a 0 .

Figure 8(a) presents a plot of the energy dependent quenching cross section computed within this 2-state model. The energy dependence of the cross section is remarkably similar to the cross sections computed with the full set of potentials. Again, we see a number of series of sharp features. The magnitude of the cross sections is larger than in the full calculations, but the size of the cross sections in the 2-state model is arbitrary and just depends upon the strength and range of the assumed coupling. within a given series are seen to increase with increasing energy. It is useful and revealing to plot the partial cross sections, namely the cross section for each partial wave of total angular momentum J for energies on and off a given feature. Figure 9 displays partial cross sections for a collision energy of 1135 cm -1 , at a sharp feature, and a slightly lower collision energy of 1110 cm -1 . For both collision energies, the partial cross sections drop abruptly to zero at the highest partial wave. Above this value of the total angular momentum J, the centrifugal barrier becomes greater than the collision energy, and the crossing region can no longer be accessed. Tunneling is not important for these systems since the reduced mass is relatively large (µ = 13.42 amu).

We see that for the sharp feature at a collision energy of 1135 cm -1 the partial cross section for the highest partial wave (J = 106) makes the dominant contribution to the integral cross section and yields a much higher integral cross section than at the slightly lower collision energy. We have examined the partial cross sections for collision energies corresponding to many of the sharp features. In all cases, the partial cross sections for the highest partial wave makes the dominant contribution to the integral cross section. The largest total angular momenta J, that make the dominant contribution, are indicated for the peaks in Fig. 8.

We can show that the sharp features are due to resonances by examining the scattering wave functions. Figure 10 displays the 1 ∆ component of the scattering wave function for several partial waves at a collision energy of 1135 cm -1 , which is the energy of one of the sharp features in Fig. 8. The bottom panel of Fig. 10 displays the potential energy curves of the 1 Σ + and 3 Π states, as well as the adiabatic potentials for total angular momentum J = 106, which makes the dominant contribution to the integral cross section at this energy. The collision energy is indicated in this panel with a horizontal line. It can be seen that this energy is very slightly smaller than the height of the centrifugal barrier.

It is interesting to compare the amplitude of the 1 ∆ component of the scattering wave function inside and outside the centrifugal barrier. For the highest, and dominant, partial wave (J = 106) that contributes to the integral cross section, the amplitude is considerably enhanced inside the barrier. This implies that the collision complex is quasi-bound, and hence there is a greater probability for quenching to the 3 P state. For a lower partial wave, J = 101, the centrifugal barrier is significantly less than the collision energy, and the amplitude of the wave (J = 106), is seen in Fig. 10 to have a very small amplitude inside the barrier. Hence, this partial wave contributes negligibly to the electronic quenching. As noted previously, the reduced mass of the OKr collision pair is large enough that tunneling is not significant.

We have examined the scattering wave functions for collision energies corresponding to other peaks in the integral cross sections. In all cases, the amplitude of the 1 ∆ component of the scattering wave function at these energies is enhanced relative to its amplitude outside the barrier. We can conclude that the peaks in the energy dependent quenching cross sections are due to resonances. The lifetimes of the quasi-bound states are proportional to the inverse of the widths of the peaks; the increasing widths of the resonances within a given series implies that the lifetime of the quasi-bound complex is correspondingly decreasing. The number of nodes of the the 1 ∆ component inside the barrier is the same (and is equal to 13) for the peaks marked in Fig. 8(b). Hence, these resonances can be described as representing different rotational levels of a quasi-bound vibrational level. We have not carried out an analysis of all the resonances in the energy-dependent In their experiment, O( 1 D) atoms were prepared by 193 nm photolysis of N 2 O, and the state distributions were determined by extrapolation to zero delay between the photolysis and probe lasers.

For quenching by Kr, they estimated the nascent state branching ratios to be j = 0/ j = 2 = 0.33 ± 0.05 and j = 1/ j = 2 equal 0.45 ± 0.07 at an estimated nascent average collision energy of 14.7 kcal mol -1 , or 5141 cm -1 . For this collision energy, we compute the following branching ratios: 1.126 for j = 0/ j = 2 and 0.073 for j = 1/ j = 2. There is hence significant disagreement between the experimentally measured and computed state branching ratios. We have no explanation for this disagreement. On the experimental side, perhaps there was more collisional relaxation of the nascent state distribution than anticipated. Also, the relative detection efficiency for each of the fine-structure levels must be determined; this was carried by interrogating a thermalized sample of 3 P j levels.

Rate Constants

We have employed the collision energy dependent cross sections for quenching of O( 1 D) to the sum over all the O( 3 P j ) fine-structure levels to determine the rate constants for O-Rg, where Rg = Ar, Kr, and Xe. The relationship between the quenching cross sections σ O( 1 D)-Rg and rate constants

k O( 1 D)-Rg is 36 k O( 1 D)-Rg (T ) = 8 πµ(k B T ) 3 1/2 ∞ 0 E c σ O( 1 D)-Rg (E c )e -E c /k B T dE c (13)
where E c is the collision energy and k B is the Boltzmann constant. The upper range of integration was taken to be E c = 5000 cm -1 .

Figure 3 presents a comparison of our computed O( 1 D)-Kr quenching rate constant with the experimental measurements of this rate constant, including the present measurements reported in Sec. II B. Our computed rate constant is significantly lower than the experimental measurements.

However, the calculation does satisfactorily describe the temperature dependence of the experimental data.

We compare in Fig. 11 Both experiment and theory display the trend of increasing magnitude of the quenching rate constants as one goes down the periodic table from Ar, to Kr, and then Xe. Two general trends were found in the potentials and spin-orbit matrix elements. Firstly, the depth of the well for the 1 Σ + state, which derives from the O( 1 D) + Rg asymptotic limit increases significantly as we go down the periodic table. More importantly, the magnitude of the spin-orbit matrix elements, in particular the coupling between the 1 Σ + state and the Ω = 0 + component of the 3 Π state, at the crossing radius increases substantially as we go down the periodic table. Similarly, the rate constant for quenching of S( 1 D) by Ar (Ref. 8) is approximately an order of magnitude larger than the rate constant for O( 1 D) quenching by Ar. This again reflects the larger size of the spin-orbit coupling in the heavier S atom.

We do find an interesting structure to the computed quenching cross sections, which is unfortunately currently not amenable to experimental observation. The energy-dependent quenching cross sections for all three collision pairs, especially O( 1 D)-Kr and O( 1 D)-Ar, show many sharp resonance peaks due to quasi-bound levels trapped by the centrifugal barrier for the highest partial wave for which quenching occurs. These resonances, which occur for only one partial wave at a given collision energy, are a consequence of the relatively large reduced mass of these collision pairs.

V. SUPPLEMENTARY MATERIAL

See the supplementary material for the file PES.tar, which contains tables of the O( 1 D, 3 P)-Ar, Kr interaction energies and spin-orbit matrix elements of the internuclear separation R.

VI. ACKNOWLEDGMENTS

The theoretical work was supported by the National Science Foundation under grant no. CHE-1565872.

2 - 5 )

 25 , Kr (Refs. 2-4), and Xe (Refs. 2 and 3) for temperatures ranging from ∼ 400 K down to 50 K. Garofalo et al. 6 investigated the dynamics of the quenching of O( 1 D) by Xe atoms in a molecular beam experiment through measurement of the differential cross section of the quenched O( 3 P) atoms. They found marked large-scale oscillations in the angle-resolved scattering of the O( 3 P) and Xe products. To interpret these intriguing observations, Dagdigian, Kłos, and Alexander 7 carried out time-independent quantum close-coupling calculations describing this quenching process, employing computed potential energy curves and spin-orbit coupling matrix elements. Comparison of the experimental and theoretical results showed that the oscillations are caused by Stueckelberg interferences, in which the quantum phases of the multiple quenching pathways interfere constructively and destructively.

  and results on the measurement of the O( 1 D) quenching rate constant by Kr. Sections III A and III B describe the ab initio calculation of the O( 1 D, 3 P)-Kr, Ar potential energy curves and spinorbit matrix elements and the details about the scattering calculations, respectively. Section III C 1 presents and interprets the computed cross sections for electronic quenching of O( 1 D) by Kr and Ar. The computed O( 3 P) fine-structure relative populations formed in the quenching by Kr are compared with a measurement by Matsumi et al. 11 in Sec. III C 2. Section III C 3 presents computed temperature-dependent rate constant for the quenching of O( 1 D) by Ar, Kr, and Xe and compares these with the experimental measurement. A discussion follows in Sec. IV.

Electronically excited state O( 1 D 2 +

 12 ) atoms [in addition to coproducts O 2 (a 1 ∆ g )] were created directly in the supersonic flow by the 10 Hz pulsed laser photolysis of ozone (O 3 ) at 266 nm with energies around 25 mJ, in a 6 mm diam beam. This process is known to occur with a quantum yield of approximately 0.9 (Ref. 15) at this wavelength, with O( 3 P) and O 2 (X 3 Σ - g ) as minor photolysis products. Ozone was generated by the ultraviolet (UV) irradiation by a high-pressure mercury lamp of pure O 2 flowed continuously into a quartz cell upstream of the reactor. The cell itself was maintained at high pressure (700 Torr), thereby promoting O 3 formation through the O( 3 P) + O M recombination reaction. The effluent of the cell (essentially O 2 with a small percentage of O 3 ) was then mixed with the carrier gas Ar and the co-reagent gas Kr and carried into the Laval nozzle reservoir. O( 1 D) atoms were probed directly in this work by pulsed vacuum ultraviolet laser induced fluorescence (VUV LIF) through the 3s 1 D -2p 1 D transition at 115.215 nm. The procedure

FIG. 1 .

 1 FIG. 1. O( 1 D) VUV LIF signal as a function of delay time recorded at 127 K. (red solid circles) without added Kr ([Ar] = 1.26×10 17 cm -3 ); (blue solid squares) with [Kr] = 1.48×10 16 cm -3 . The solid red and blue lines represent the best fit derived by fitting a single exponential decay function to the experimental data points. Pseudo-first-order rates were extracted from the time constant of the fit.

1 FIG. 3 .

 13 FIG. 3. Rate constants for the O( 1 D) + Kr quenching reaction as a function of temperature. Experimental results: red filled circles: this work; blue squares, Davidson et al. (Ref. 2); green triangle, Shi and Barker (Ref. 3); cyan diamonds, Blitz et al. (Ref. 4). Computed value: black line, this work.

Figure 4

 4 Figure 4 displays the potential energy curves for all the O( 3 P, 1 D)-Kr, Ar states. As in the case of O-Xe, 7 we see that the 1 Σ + state is the most strongly bound singlet state for both complexes. For this state, the O 2p orbital aligned along the internuclear axis is vacant. The binding energy of this state is the smallest for O-Ar, significantly larger for O-Kr, and larger still for O-Xe [D e = 13365.8 cm -1 ]. 7 The other electronic states have shallow van der Waals wells. The 1 ∆ and 3 Σ -states are the most repulsive; for these states the O 2p orbital aligned along the internuclear axis is doubly occupied. The values of the dissociation energies D e and equilibrium internuclear separations R e for the computed O-Kr and O-Ar triplet and singlet states are listed in TableII. We compare the val-

FIG. 5 .

 5 FIG. 5. Spin-orbit matrix elements between OKr (top panels) and OAr (bottom panels) wave functions. Matrix elements between 3 P projection states: A xy ≡ 3 Π x |L z S z | 3 Π y and A zy ≡ 3 Π x |L x S x | 3 Σ .

Figure 6

 6 Figure 6 displays the integral cross sections for formation of the individual 3 P j fine-structure levels, as well as the sum of these fine-structure dependent cross sections, as a function of collision energy for the quenching of O( 1 D) by Kr and Ar. The cross sections show slow oscillations in magnitude as a function of the collision energy, as well as superimposed sharp features. This behavior is similar to what we observed for the energy dependence of the O-Xe quenching cross sections, 7 except that the sharp features are much more prominent for O-Kr and O-Xe. It can be seen that over the range of collision energies investigated only the 3 P j , j = 2 and 0, fine-structure states are significantly populated in the quenching of O( 1 D) by Kr and Xe. This is in contrast to the quenching of O( 1 D) by Xe,7 for which the 3 P 1 state is formed with significant yield for collision energies greater than approximately 2000 cm -1 . For this system, quenching to form the 3 P 1 state can occur through the crossing of the 1 Π and 3 Σ -states. As discussed in Sec. III A, this crossing lies too high to be accessible in O-Kr and O-Ar collisions.

2 FIG. 6 .

 26 FIG. 6. Integral cross sections for the quenching of O( 1 D) by (a) Kr and (b) Ar to the individual O( 3 P j )fine-structure as a function of the collision energy. Also plotted is the cross section summed over the finestructure levels.

Figure 8 ( 2 FIG. 7 .

 827 Figure 8(b) displays an expanded plot over the narrow energy range. The widths of the features

FIG. 8 .

 8 FIG. 8. (a) Integral cross sections for the quenching of O( 1 D) by Kr as a function of the collision energy within the 2-state model. (b) Plot with an expanded collision energy range; the numbers above the peaks indicate the total angular momentum J of the dominant, and highest, partial wave contributing to the integral cross section at that energy.

1 ∆ 1 FIG. 9 .

 119 FIG. 9. (a) Partial cross sections as a function of the total angular momentum J for the quenching of O( 1 D) by Kr as a function of the collision energy within the 2-state model at collision energies of (a) 1135 cm -1 and (b) 1110 cm -1 .

2 FIG. 10 .

 210 FIG.10. Top three panels: The 1 ∆ component of the scattering wave function for several partial waves (J = 101, 106, and 107) for a collision energy of 1135 cm -1 . Bottom panel: Potential energy curves for the 1 Σ + and 3 Π states in the 2-state model and adiabatic energies E 1 and E 2 for total angular momentum J = 106. Also plotted in the bottom panel is the total energy for a collision energy E c = 1135 cm -1 .

2 .

 2 Comparison with VUV Experiment Matsumi et al. 11 utilized VUV LIF in a cell to determine the O( 3 P j ) fine-structure state distribution in the quenching of O( 1 D) in collisions with several atomic and molecular collision partners.

FIG. 11 .

 11 FIG. 11. Experimentally measured and computed rate constants for the quenching of O( 1 D) by (a) Ar and (b) Xe. Experimental values in panel (a): Red squares, Grondin et al. (Ref.5); blue squares, Davidson et al. (Ref. 2); green triangle, Shi and Barker, upper limit (Ref. 3); cyan diamonds, Blitz et al. (Ref. 4). Experimental values in panel (b): Blue squares, Davidson et al. (Ref. 2); green triangle, Shi and Barker (Ref. 3). Computed values: black lines, this work.

  Measured pseudo-first-order rate constants as a function of excess Kr concentration. Red open circles: data recorded at 296 K; blue solid squares: data recorded at 50 K. The solid red and blue lines represent weighted fits to the respective data sets (see text for details), yielding the second-order rate constant from the slope.

	The O( 1 D) signal decays rapidly even in the absence of Kr [k ′ = (7.58 ± 0.19) × 10 4 s -1 for the
	red data points shown in Fig. 1] due to the efficient relaxation of O( 1 D) atoms by the carrier gas
	Ar with a roughly temperature independent rate constant k O( 1 D)-Ar = 6.0 × 10 -13 cm 3 s -1 (Refs. 2
	and 5). In this particular example, [Ar] = 1.26 × 10 17 cm -3 [see Table 1 of Grondin et al. (Ref. 5)],
	so that k O( 1 D)+Ar [Ar] = 7.56 × 10 4 s -1 . The contributions of processes (4) -(6) are clearly quite
	negligible with respect to processes (2) and (3) in this example and in all other decays, as already
	argued by Grondin et al.

Some typical O( 1 D) decay curves recorded at 127 K in the presence and absence of Kr are shown in Fig. 2.

FIG. 2.

TABLE II .

 II Computed dissociation energies [D e , in cm-1 ] and equilibrium internuclear separations [R e , in a 0 ] for the electronic states of OKr and OAr.

	State		OKr		OAr	Empirical OKr a	Empirical OAr a
		D e	R e	D e	R e	D e	R e	D e	R e
						O( 3 P) + Rg
	3 Σ -	54.3	7.453	46.3	7.201	71.0 ± 7.1	7.65 ± 0.15 58.1 ± 5.8	7.28 ± 0.14
	3 Π	100.2	6.540	82.7	6.349	105.7 ± 10.6 6.75 ± 0.14 83.9 ± 8.4	6.52 ± 0.13
						O( 1 D) + Rg
	1 Σ +	7150.7 3.476	3129.2 3.366	-	-	-	-
	1 Π	81.7	6.615	61.6	6.524	-	-	-	-
	1 ∆	54.4	7.342	46.9	7.082	-	-	-	-

a Potentials derived from the scattering experiments of Aquilanti et al. (Ref. 23).


Ousmane Diop
email: ousmane.diop@inrae.fr

Madalena Chaves

Laurent Tournier

Qualitative analysis of mammalian circadian oscillations: cycle dynamics and robustness

Keywords: Mammalian circadian clock, Asynchronous Boolean network, Complex attractor, Summary graph

In asynchronous Boolean models, periodic solutions are represented by terminal strongly connected graphs, which are typically composed of hundreds of states and transitions. For biological systems, it becomes a challenging task to compare such mathematical objects with biological knowledge, or interpret the transitions inside an attractor in terms of the sequence of events in a biological cycle. A recent methodology generates summary graphs to help visualizing complex asynchronous attractors and order the dynamic progression based on known biological data. In this article we apply this method to a Boolean model of the mammalian circadian clock, for which the summary graph recovers the main phases of the cycle, in the expected order. It also provides a detailed view of the attractor, suggesting improvements in the design of the model's logical rules and highlighting groups of transitions that are essential for the attractor's robustness.

Introduction

The analysis of periodic orbits and their properties remains a most challenging problem in dynamical systems theory. Many living systems exhibit periodical dynamics and the current literature covers a large diversity of mathematical models used to represent, explore, and study the mechanisms leading to physical or biological rhythms [START_REF] Forger | Biological clocks, rhythms, and oscillations: the theory of biological timekeeping[END_REF]. A thorough analysis of such cyclic attractors opens the door to a whole family of meaningful questions related to the robustness of the oscillatory behavior, the estimation and control of the period or amplitude of oscillations in terms of the parameters of the system, the location of the orbit in the state space, etc.

Very little is known on how to express the properties of a periodic solution in terms of the system's parameters but qualitative models, such as piecewise linear or Boolean models, suggest some ideas. Piecewise linear systems partition the state space into regions where solutions of the system can be explicitly computed leading, in some examples, to the estimation of the period and other quantities in terms of the parameters [START_REF] Ndiaye | Oscillations induced by different timescales in signal transduction modules regulated by slowly evolving protein-protein interactions[END_REF][START_REF] Poignard | A stability result for periodic solutions of nonmonotonic smooth negative feedback systems[END_REF]. Boolean models provide an ideal framework to analyze qualitative dynamical properties by enabling algorithmic approaches to characterize, for instance, the location of a periodic orbit in the state space, or the influence of the interaction graph on the dynamical behavior [START_REF] Tournier | Uncovering operational interactions in genetic networks using asynchronous Boolean dynamics[END_REF][START_REF] Remy | Boolean dynamics of compound regulatory circuits[END_REF].

In an asynchronous Boolean network, a periodic orbit corresponds to a terminal strongly connected component of the state transition graph. If the computation of such an object is not an issue form a theoretical point of view, its size can grow very large, strongly limiting the biological interpretation of the states and transitions inside the attractor. To tackle this issue, general approaches can be used such as model checking techniques [START_REF] Traynard | Model-based investigation of the circadian clock and cell cycle coupling in mouse embryonic fibroblasts: Prediction of reverb-α up-regulation during mitosis[END_REF]. Recently, a dedicated method developed by Diop et al. [START_REF] Diop | Summarizing complex asynchronous Boolean attractors, application to the analysis of a mammalian cell cycle model[END_REF] proposes to generate a summary graph of an asynchronous attractor, based on a classification of its states according to experimentally observed phases of the biological system. The summary graph provides a qualitative view of the general progression along the periodic orbit, capturing the underlying dynamics within the cyclic attractor. In [START_REF] Diop | Summarizing complex asynchronous Boolean attractors, application to the analysis of a mammalian cell cycle model[END_REF] it was successfully applied to a Boolean model of the mammalian cell cycle [START_REF] Fauré | Dynamical analysis of a generic boolean model for the control of the mammalian cell cycle[END_REF].

In this paper, we propose to apply the summary graph method to the mammalian circadian clock, a biological rhythm which is based on the interactions among five main proteins CLOCK:BMAL1, REV-ERBα, ROR, PER, and CRY. The core of the clock mechanism is formed by three feedback loops. First, the BMAL1 complex promotes the transcription of Per and Cry genes. The corresponding proteins bind to form a complex PER:CRY, which then translocates to the nucleus where it will block the transcriptional activity of BMAL1, thus forming a first negative feedback loop. In addition, BMAL1 also promotes the transcription of the two genes Rev-erb and Ror. Eventually, Bmal1 transcription will be inhibited by REV-ERB and activated by ROR, leading to a second negative loop and a new positive loop, respectively. The periodic behavior of the clock system is determined in large part by the phase opposition between CLOCK:BMAL1 and PER:CRY, which corresponds to the day/night succession (BMAL1 peaks during the day).

Examples of Boolean models for circadian clocks include one for the plant Arabidopsis thaliana [START_REF] Akman | Digital clocks: simple boolean models can quantitatively describe circadian systems[END_REF] and a compact mammalian clock model [START_REF] Comet | Simplified models for the mammalian circadian clock[END_REF] which reproduces the interplay between the negative feedback loops induced by BMAL1 activity. However, our present objective is to have a deeper understanding of the circadian cycle recently developed by Almeida et al. [START_REF] Almeida | Transcription-based circadian mechanism controls the duration of molecular clock states in response to signaling inputs[END_REF]. This is a continuous model of the core clock mechanism that faithfully reproduces the circadian rhythm, by including not only the five main proteins but also their transcription regulated by clock controlled elements (CCE). To use the graph method [START_REF] Diop | Summarizing complex asynchronous Boolean attractors, application to the analysis of a mammalian cell cycle model[END_REF], we first construct a Boolean version of the continuous model in [START_REF] Almeida | Transcription-based circadian mechanism controls the duration of molecular clock states in response to signaling inputs[END_REF] that exhibits one cyclic attractor and, in addition, correctly reproduces the effect of some well known gene knock-outs. Next, Section 3 analyzes the cyclic attractor with the summary graph method, by classifying groups of Boolean states according to their corresponding circadian time zones. Finally, in Section 4 the summary graph is further used to identify key groups of transitions within the attractor and relate them back to parts of the network's logical rules. This analysis has two outcomes, first suggesting an improved Boolean rule for one of the variables where the model lacked clear information. As a second outcome, our analysis predicts that some transitions, while allowing some short-circuits between phases, seem to contribute to the overall robustness of the attractor and globally ensure the good progression of the cell clock.

Proposing a new Boolean model of the circadian clock

As evoked previously, the following Boolean model of the mammalian circadian clock is highly inspired by the continuous model developed in [START_REF] Almeida | Transcription-based circadian mechanism controls the duration of molecular clock states in response to signaling inputs[END_REF]. In addition to the three feedback loops already described, this model takes into account the transcription of the five clock proteins, each regulated through a particular combination of transcription factors. The latter bind to specific sites on the promoters, called clock controlled elements (CCEs): Ebox (enhancer box), Dbox (DBP/E4BP4 response element), and Rbox (REV-ERB/ROR response element). With the introduction of Dbox (which activates both REV-ERB and PER), two new proteins are added to the model, DBP and E4BP4, each also regulated by one of the CCE.

A complete and detailed justification of the continuous model assumptions and construction can be found in [START_REF] Almeida | Transcription-based circadian mechanism controls the duration of molecular clock states in response to signaling inputs[END_REF], but we provide a brief summary of the eight variables and corresponding differential equations in Table 1.

Construction of the Boolean model

The design of a Boolean or of a continuous model for the same biological system differs in some fundamental aspects and, in general, there is no direct equivalence between terms in the two frameworks. For instance, activation and inhibition links typically have clear logical representations, but the effects of detailed mass-action kinetics or mass conservation laws are harder to represent in a Boolean model, and may require the definition of new variables. The purpose of our Boolean model is to transcribe as closely as possible the interactions in the continuous model in [START_REF] Almeida | Transcription-based circadian mechanism controls the duration of molecular clock states in response to signaling inputs[END_REF], as described below. The continuous equations and corresponding logical rules are shown in Table 1, for a clear comparison between the two models.

In [START_REF] Almeida | Transcription-based circadian mechanism controls the duration of molecular clock states in response to signaling inputs[END_REF], there are three CCE named Dbox, Ebox, and Rbox, each responding to the conjugation of two components, an activator and an inhibitor. More precisely, Dbox is activated by DBP and inhibited by E4, Ebox is activated by Bmal and inhibited by CRY and Rbox is activated by ROR and inhibited by REV. The continuous equations are formed by a synthesis term depending on one or two CCE and a degradation term. The CRY , P ER and P C equations also contain the binding and dissociation terms denoted M P C in Table 1.

In general, Boolean variables are assumed to degrade when not updated in the next step and, therefore, linear degradation terms do not appear explicitly. Thus, Table 1. Differential equations of the circadian clock model in [START_REF] Almeida | Transcription-based circadian mechanism controls the duration of molecular clock states in response to signaling inputs[END_REF] and corresponding logical rules in the Boolean model. The term MP C = γpcP ER•CRY -γcpP C represents the kinetics of the complex binding: PER+CRY PC (mass action law). We use classical Boolean operators: X ∨ Y (X or Y ), X ∧ Y (X and Y ) and X (not X).

Continuous equation

Logical rule

dBM AL1/dt = Rbox -γ bp BM AL1 • P C Bmal = Rbox ∧ P C dROR/dt = Ebox + Rbox -γrorROR ROR = Rbox2 ROR2 = Ebox ∧ Rbox ∧ ROR dREV /dt = 2Ebox + Dbox -γrevREV REV = Ebox ∨ REV 2 REV 2 = Ebox ∧ Dbox ∧ REV dDBP/dt = Ebox -γ dbp DBP DBP = Ebox dE4/dt = 2Rbox -γe4E4 E4 = Rbox dCRY /dt = Ebox + 2Rbox -MP C -γcCRY CRY = Ebox ∨ CRY 2 CRY 2 = Ebox ∧ Rbox ∧ CRY dP ER/dt = Ebox + Dbox -MP C -γpP ER P ER = Ebox ∨ Dbox dP C/dt = MP C -γ bp BM AL1 • P C P C = P ER ∧ CRY
a naive approach to construct a corresponding Boolean model is to combine the synthesis terms as logical conjugations or disjunctions of the given variables and directly obtain the rule, for instance: DBP = Ebox or E4 = Rbox. However, not all continuous equations follow this simple construction and other properties that strongly contribute to the dynamics must be taken into account. Indeed, (i) some variables have nonlinear degradation (cf. BM AL1), (ii) others contain mass-action terms (cf. CRY ,P ER,P C) and (iii) three of the variables are regulated by two CCEs and are themselves regulators of CCEs (ROR, REV , CRY ). This last property implies the existence of different thresholds for the different regulatory activities. For this reason, ROR, REV , and CRY are assumed to have an extra discrete level, here represented by extra Boolean variables ROR2, REV 2, and CRY 2 (following [START_REF] Van Ham | How to deal with variables with more than two levels[END_REF]), under the assumption

x =    0 when X = X2 = 0, 1 when X = 1, X2 = 0, 2 when X = X2 = 1.
The states corresponding to X = 0, X2 = 1 have no biological meaning (also called "forbidden") and trajectories towards these states from the other "biological" states must be excluded. This can be achieved by using the method described in [START_REF] Chaves | Comparing Boolean and piecewise affine differential models for genetic networks[END_REF] which complements the rules for

X = (• • •) and X2 = (• • •) as follows: X = (• • •) ∨ X2 and X2 = (• • •) ∧ X.
Since the CCEs appear in additive form in the continuous equations, we assume that only one CCE is sufficient to trigger the activity of X, while the two CCEs are needed to trigger the activity of X2. The order in which the CCEs are activated and in turn activate each variable, was decided by comparing to the continuous solutions. The expressions for the CCEs are written as logical conjunctions:

Ebox = Bmal ∧ CRY 2, Dbox = DBP ∧ E4, Rbox = ROR ∧ REV , Rbox 2 = ROR2 ∧ REV .
In the continuous model, there is only one nonlinear degradation term, corresponding to the inhibition of BMAL1 transcription by PC (-γ bp BM AL1 • P C). Accordingly, we assume Bmal is explicitly repressed by P C but, conversely, PC is not strongly affected (indeed, our analysis showed that a rule of the form P C = P ER ∧ CRY ∧ Bmal prevents oscillatory behavior).

The binding of P ER and CRY to form the complex P C is described by mass-action kinetics in the continuous model (see term M P C in Table 1). In the Boolean model, for simplicity, we assumed that P C is produced when both P ER and CRY are available, leading to the rule P C = P ER ∧ CRY , but no explicit effect from P C on P ER or CRY . Although there is no systematic way to translate mass-action kinetics or other mass conservation laws into Boolean factors, and these are usually treated on a case-by-case basis, we will see in Part 4 that the analysis by the summary graph method suggests a refinement of the P ER rule. This refinement can be interpreted as a more suitable way to include the mass-action terms into the circadian clock Boolean model.

In this way, we obtained a qualitative multi-valued model which closely translates the differential model of [START_REF] Almeida | Transcription-based circadian mechanism controls the duration of molecular clock states in response to signaling inputs[END_REF]. The multi-valued model is equivalent to the 11-dimensional Boolean network depicted in Table 2, to which the methodology in [START_REF] Diop | Summarizing complex asynchronous Boolean attractors, application to the analysis of a mammalian cell cycle model[END_REF] can now be applied.

Bmal = ROR ∧ REV ∧ P C ROR = ROR2 ∨ REV REV = (Bmal ∧ CRY 2) ∨ REV 2 DBP = Bmal ∧ CRY 2 E4 = ROR ∧ REV CRY = Bmal ∨ CRY 2 P ER = (Bmal ∧ CRY 2) ∨ (DBP ∧ E4) P C = P ER ∧ CRY REV 2 = Bmal ∧ CRY 2 ∧ DBP ∧ E4 ∧ REV CRY 2 = Bmal ∧ CRY 2 ∧ ROR ∧ REV ∧ CRY ROR2 = Bmal ∧ CRY 2 ∧ ROR ∧ REV

First dynamical analysis of the Boolean model

To compute and analyze the dynamical behavior of a Boolean model X i = F i (X 1 , . . . , X n ), i = 1, . . . , n, we need to specify an updating order for the variables. Applying the logical rules in the order defined by the updating schedule leads to a state transition graph (STG) with 2 n states, where a sequence of transitions represents a trajectory of the system. To analyze the STG, we first compute its strongly connected components (SCCs) which are defined as sets of states C, such that for every pair of states x, y ∈ C, there exist two paths (or sequences of transitions) in C leading from x to y and from y to x. SCCs may consist of single or multiple states and may have incoming and outgoing transitions, but two distinct SCCs can not be mutually connected, otherwise they would form a single SCC. The asymptotic behavior of the system is thus characterized by the SCCs without outgoing transitions, also called terminal SCCs or attractors. An attractor with multiple states represents a periodic orbit of the system.

As a preliminary analysis, we considered the basic synchronous updating schedule, where all variables are simultaneously updated:

X i [t + 1] = F i (X[t]).
The synchronous STG of the model contains a single attractor, composed of only five states. This simple cyclic attractor captures the Bmal/PER:CRY phase opposition, a central feature in the circadian clock (see Figure 1). In this paper, we will prefer an asynchronous updating schedule as it is much more realistic from a biological perspective. In this scheme, at most one variable is updated at each instant. To construct the STG for this asynchronous schedule, for each state X = (X 1 , . . . , X n ) define the subset of variables I X = {i : X i = F i (X)}. Then, for each i ∈ I X , add a transition X → Y where Y i = F i (X) and Y j = X j for all j = i. More details on the asynchronous strategy can be found for instance in [START_REF] Tournier | Uncovering operational interactions in genetic networks using asynchronous Boolean dynamics[END_REF]. From a biological point of view, asynchronous updating is preferable as it allows for variability and different timescales in the network interactions. From a graph theoretical point of view, synchronous and asynchronous STGs are very different as the asynchronous STG is non deterministic, i.e. a state no longer has one and only one successor. Usually, this generates lots of transitions, making biological interpretation more difficult.

When applied to our circadian model, the asynchronous strategy still gives a single attractor, but the terminal SCC is now composed of 442 states and 1737 transitions. Interestingly, the five states of the synchronous attractor are all present in the asynchronous attractor, but the simple cyclic trajectory has now been replaced by a complex graph with a much larger amount of details. To further analyze this attractor and validate our model, we consider six different versions of the model representing six known mutations (see Table 3). For some of them, circadian oscillations may be completely lost, which in our model translates to the attractor shrinking to a single fixed state. For others, oscillations endure but are somewhat degraded (eg. with a shorter period); in those cases our model conserves a complex attractor, but with less states than the original. These results seem to confirm that the model in Table 2 reproduces the essential core of the circadian clock dynamics.

Table 3. Effects of some mutations on the dynamics of the Boolean model.

Mutation Biological phenotype

Effect on the attractor Bmal=0 Arrhythmic (complete loss of circadian oscillations) [START_REF] Ko | Molecular components of the mammalian circadian clock[END_REF][START_REF] Ripperger | The daily rhythm of mice[END_REF][START_REF] Rosensweig | Periodicity, repression, and the molecular architecture of the mammalian circadian clock[END_REF] Single state attractor 0100100000 PER=0 Abnormal circadian oscillations [START_REF] Ko | Molecular components of the mammalian circadian clock[END_REF][START_REF] Ripperger | The daily rhythm of mice[END_REF][START_REF] Rosensweig | Periodicity, repression, and the molecular architecture of the mammalian circadian clock[END_REF] Attractor of 114 states CRY=0 Abnormal circadian oscillations [START_REF] Ko | Molecular components of the mammalian circadian clock[END_REF][START_REF] Ripperger | The daily rhythm of mice[END_REF][START_REF] Rosensweig | Periodicity, repression, and the molecular architecture of the mammalian circadian clock[END_REF] Attractor of 106 states REV=0 Shorter cycle period [START_REF] Preitner | The orphan nuclear receptor rev-erbα controls circadian transcription within the positive limb of the mammalian circadian oscillator[END_REF][START_REF] Ripperger | The daily rhythm of mice[END_REF] Attractor of 80 states ROR=0 Arrhythmic [START_REF] Ko | Molecular components of the mammalian circadian clock[END_REF][START_REF] Ripperger | The daily rhythm of mice[END_REF] Single state attractor 00000000000 REV=1 Arrhythmic [START_REF] Kornmann | System-driven and oscillator-dependent circadian transcription in mice with a conditionally active liver clock[END_REF] Single state attractor 00100000000

Mutant analysis constitutes an interesting way to validate a discrete dynamical model. Indeed, the comparison between wild type and mutant phenotypes usually provides qualitative differences, such as the disappearance of oscillations for instance, that can be well captured by a Boolean model. Nevertheless, when dealing with such a complex attractor (hundreds of states, thousands of transitions), a more direct comparison of the attractor with biological data is rapidly limited, hindering model validation. For example, in Table 3 the degradation of circadian oscillations is paralleled with the number of states in the attractor, which is questionable. In the following we use the methodology proposed in [START_REF] Diop | Summarizing complex asynchronous Boolean attractors, application to the analysis of a mammalian cell cycle model[END_REF], constructing a reduced version of the attractor based on biological knowledge. This summary graph leads to refine the analysis of the attractor, providing further validation of the model.

Comparing the attractor with circadian oscillations

Dividing the circadian cycle into qualitative phases

In order to further analyze the model's attractor, a necessary first step is to classify its states into groups that will correspond to different stages of the circadian clock. This is essentially a modeling step, therefore there is not a unique way to make this classification, as it is based on a compromise between available biological data on the one hand, and the different variables and interactions included in the model on the other. In the following, we give a brief description of the main regulatory events during the circadian cycle, together with the modeling choices we made to deduce the corresponding partition of the attractor's states. For a comprehensive biological review we mainly referred to [START_REF] Takahashi | Transcriptional architecture of the mammalian circadian clock[END_REF]. Note that this article is based on data at the transcriptional level; therefore, we sometimes used other sources to complete our classification (see [START_REF] Almeida | Transcription-based circadian mechanism controls the duration of molecular clock states in response to signaling inputs[END_REF][START_REF] Relógio | Tuning the mammalian circadian clock: robust synergy of two loops[END_REF] and references therein).

As already mentioned, a hallmark of circadian rhythm progression is the phase opposition between the CLOCK:BMAL1 complex on the one hand and the PER:CRY complex on the other. This opposition divides the cycle into two major steps, approximately correlated with the day/night separation (see Figure 2). More precisely, [START_REF] Takahashi | Transcriptional architecture of the mammalian circadian clock[END_REF] introduces two biological phases, respectively called Activation and Repression. The first one takes place during the day and corresponds to the activation of Bmal, while PER:CRY is absent. The second one sees the repression of Bmal and takes place during the night. Projecting this on the variables of our model, this leads to consider two groups of states: one where Bmal = 1, P C = 0 and one where Bmal = 0, P C = 1. Fig. 2. Description of the main qualitative stages of the circadian clock. CT stands for Circadian Time and is an standard marker of time arbitrarily starting (CT0) at the beginning of activity for a diurnal organism. Top: temporal succession of biological phases as described in [START_REF] Takahashi | Transcriptional architecture of the mammalian circadian clock[END_REF]; middle: projection of the main regulatory events on the model's variables; bottom: corresponding qualitative phases Vi defined in (1).

These two main phases are separated by intermediate phases called Transcription and Poised state-Derepression in [START_REF] Takahashi | Transcriptional architecture of the mammalian circadian clock[END_REF]. In the Transcription phase, Per and Cry genes are transcribed, followed by their complexation and the translocation of the complex into the nucleus. With the lack of precise timing of these events, we decided to subdivide this into a first step where PER or CRY are present but not at the same time, followed by a second step where they are both present. The latter corresponds to the pre-formation step of the PER:CRY complex, i.e. PER and CRY have sufficiently accumulated and the complex is about to form in the cytosol [START_REF] Rosensweig | Periodicity, repression, and the molecular architecture of the mammalian circadian clock[END_REF]. Finally, after the repression phase, the PER:CRY complex disappears due to auto-repression of PER and CRY [START_REF] Gallego | Post-translational modifications regulate the ticking of the circadian clock[END_REF]. Again, without knowing the precise timing of the disappearance of PER, CRY and the complex, we simply consider intermediary states where Bmal = P C = P ER = CRY = 0, just before the activation of Bmal and the subsequent beginning of a new cycle.

This description leads to the definition of five groups of states in the attractor, or qualitative phases, denoted by (V i ) 1≤i≤5 and defined as follows.

           V 1 = (Bmal ∧ P C) ∧ (P ER ∨ CRY ) (afternoon), V 2 = (Bmal ∧ P C) ∧ P ER ∧ CRY (late afternoon), V 3 = Bmal ∧ P C ∧ P ER ∧ CRY (transition day-night), V 4 = Bmal ∧ P C (night), V 5 = Bmal ∧ P C ∧ P ER ∧ CRY (late night to next morning). (1)
These phases are defined by taking into account the main variables of our Boolean model. They are based on qualitative and not temporal considerations; however, thanks to the description in [START_REF] Takahashi | Transcriptional architecture of the mammalian circadian clock[END_REF] we were able to approximately place them along the circadian time scale (see Figure 2).

Remark 1. For the sake of convenience, we use the same symbol V i to designate both the subset of states in the attractor and the Boolean formula describing those states. For instance, V 4 = Bmal ∧ P C denotes the set of states in the attractor such that Bmal = 0 and P C = 1.

Construction of the summary graph

Let A = (V, E) denote the attractor, which is a directed graph over |V | = 442 states. The sets V i defined by (1) are, by definition, subsets of V that are mutually exclusive. To complete them into a partition of V , introduce the set U = V \ 5 i=1 V i , containing "unclassified" states. The first thing to note with this partition is that every V i is actually not empty, confirming the attractor accurately captures all important phases of the circadian clock. To be more precise, |V We now briefly recall the definition of a summary graph (interested reader may refer to [START_REF] Diop | Summarizing complex asynchronous Boolean attractors, application to the analysis of a mammalian cell cycle model[END_REF] for more details).

Definition 1. The summary graph of the graph A = (V, E) on a partition P = {V 1 , V 2 , . . . , V k } of V is the directed graph G = (P, E), whose vertices are the V i , i = 1, . . . , k and where there is an edge from V i to V j iff i = j and there exist x ∈ V i and y ∈ V j such that (x, y) ∈ E.

From a graph theoretical point of view, G is simply the quotient graph of A on the partition P. For the sake of simplicity, we use the same symbol V i to design the set of states V i and the vertex V i of the summary graph. This enables an easy description of sets of trajectories in the attractor A: actually, each edge of the summary graph is unequivocally associated with a subset of asynchronous transitions in A. The summary graph of the attractor is depicted in Figure 3 (top left). The summary graph in Figure 3(top left) shows that the expected succes-

sion of phases V 1 → V 2 → V 3 → V 4 → V 5 → V 1 is
actually present, confirming the existence of accurate circadian oscillations in the attractor, with respect to regulators Bmal, PER , CRY and PER:CRY. Moreover, three one-directional transitions V 2 → V 3 , V 3 → V 4 , and V 4 → V 5 even indicate irreversible progression through the cycle. The edge V 2 → V 3 corresponds to the formation and translocation of the PC complex, following the accumulation of PER and CRY. This irreversibility is in adequacy with [START_REF] Rosensweig | Periodicity, repression, and the molecular architecture of the mammalian circadian clock[END_REF]: as the main role of P C is to inhibit Bmal, it will switch off only after Bmal = 0. Similarly, the edge from V 3 to V 4 corresponds to transitions where Bmal is repressed, which is known to be irreversible. As for the edge V 4 → V 5 , it is associated to a set of transitions where P C switches off and can be viewed as a consequence of PER and CRY auto-repression [START_REF] Gallego | Post-translational modifications regulate the ticking of the circadian clock[END_REF]. From a graph point of view, the presence of these three one-directional edges imposes a general orientation for the cycle, consistent with what is known on the opposition between the Activation and Repression phases of [START_REF] Takahashi | Transcriptional architecture of the mammalian circadian clock[END_REF].

In addition to the summary graph in wild type condition, Figure 3 also shows the summary graphs for the three mutants that exhibited oscillatory behaviors in Table 3: P ER = 0, CRY = 0 and REV = 0. Interestingly, the mutation of REV has a limited impact on the succession of the phases, showing this perturbation does not impair the cycle progression, as it was observed in [START_REF] Preitner | The orphan nuclear receptor rev-erbα controls circadian transcription within the positive limb of the mammalian circadian oscillator[END_REF][START_REF] Ripperger | The daily rhythm of mice[END_REF]. For P ER and CRY mutants, the perturbation is stronger and the cycle is impaired. However, we still observe a Bmal oscillation between V 1 and V 5 (notably U ), confirming circadian oscillations may persist, although in a degraded form [START_REF] Ko | Molecular components of the mammalian circadian clock[END_REF][START_REF] Ripperger | The daily rhythm of mice[END_REF][START_REF] Rosensweig | Periodicity, repression, and the molecular architecture of the mammalian circadian clock[END_REF].

Whether in wild type or in mutant conditions, the summary graph reveals to be a powerful tool to analyze asynchronous attractors of complex oscillatory biological systems, such as the circadian clock. In the next part, we propose an extension of this analysis, where we show how to use this tool to help redesign the network, and assess the general robustness of the model.

Advanced analysis of the attractor

The summary graphs in Figure 3 allow a direct comparison of the model's attractor with biological knowledge on the circadian cycle, both in wild type and mutant conditions. In particular, we were able to track the expected succession of qualitative phases within the attractor, thus validating that our model captures the essential core of the regulatory network. Nevertheless, the summary graphs also point to spurious transitions to (and from) the set of unclassified states U , responsible for "short-circuits" between phases. Such spurious trajectories may reveal abnormal behaviors, with unwanted transitions in the model, or they may catch important links ensuring the overall robustness of the cycle. In both cases, the summary graph provides an ideal tool to detect those spurious transitions and analyze them in a biological context.

Adjustment of the attractor and refinement of the model

We start with a more in-depth examination of the set U . In part 3 it is defined as 11 and is a subset of 137 states in the attractor that we were not able to classify according to the general description of the circadian cycle. By assigning the value 1 to the states in U and the value 0 to the states in V \U , we obtain a partial Boolean function (PBF) that we can identify using Boolean inference techniques (see [START_REF] Tournier | Uncovering operational interactions in genetic networks using asynchronous Boolean dynamics[END_REF][START_REF] Tournier | Optimal resource allocation enables mathematical exploration of microbial metabolic configurations[END_REF]). The inference consists in identifying minimal supports 3 , i.e. minimal subsets of variables that are sufficient to reproduce the PBF. This function provides an easy way to characterize the unclassified states, showing they can be decomposed into two components:

U := V \ ( i V i ) ⊆ {0, 1}
U = Bmal ∧ P C ∧ (P ER ∨ CRY ) U1 ∨ Bmal ∧ P C ∧ P ER U2 . (2)
While the inference of a PBF does not lead to a unique solution in general, here (2) is actually unique, as the set {Bmal, PC, PER, CRY} is the only minimal support of the PBF U (it will be the case for all PBF identified in this article). As before, we conveniently use the same symbol U to designate the set of unclassified states (within the attractor) and the Boolean function that characterizes them (see Remark 1 above). The Boolean formula of U is decomposed as in (2) to highlight two separate subsets of unclassified states: the first one U 1 is composed of states where variables Bmal and P C are both off whereas the second one U 2 is composed of states where variables Bmal and P C are both on. Specifically, in this second subset the PER:CRY complex is forming (since it has not repressed Bmal yet) while PER is switched off. Clearly, such states should not exist in the attractor since PER is essential for complex formation. In order to observe the interplay between U 1 , U 2 with the different phases V i , we reconstruct the summary graph of the attractor, using Definition 1 on the new partition {V 1 . . . , V 5 , U 1 , U 2 } (Figure 4, left).

By looking at this graph, one can see that the set U 2 is connected to the rest of the attractor mainly through phase V 3 . More precisely, only six transitions in the attractor are responsible for the entry into U 2 , all coming from V 3 . Since an asynchronous strategy is used, all these transitions can be traced back to situations where PER has disappeared before the complex PER:CRY has completed the repression of Bmal, leading to a contradiction. Therefore, we decide to suppress these transitions in the attractor, leading to the disconnection of the set U 2 (Figure 4, center). After the removal of the six transitions in the graph of the attractor, we obtain a new amended model by applying a few steps. First, we recompute the strongly connected component decomposition of the truncated attractor. We find a unique terminal SCC of 393 states, in which the states in U 2 have disappeared.

From this SCC we reconstruct the partial truth table of the network on these 393 states, and then use an inference technique to identify the logical rules corresponding to this partial truth table. The last two steps (construction of partial truth table, followed by inference) are described in more details in [START_REF] Diop | Summarizing complex asynchronous Boolean attractors, application to the analysis of a mammalian cell cycle model[END_REF] and, in a different context, in [START_REF] Tournier | Uncovering operational interactions in genetic networks using asynchronous Boolean dynamics[END_REF]. Finally, the only affected rule of the network is the rule of variable P ER, with the addition of a new (unique) clause:

P ER = (Bmal ∧ CRY 2) ∨ (DBP ∧ E4)∨(Bmal ∧ PC). (3)
To verify that this modification did not alter other parts of the dynamics, we made the same analyses as in Parts 2 and 3, confirming that the slight modification in (3) was sufficient to get rid of unwanted states U 2 while conserving good dynamical properties, in wild type and mutant conditions. The summary graph of this new amended model is depicted in Figure 4, on the right. The modified rule (3) adds a single new interaction in the network, which is a positive effect of P C onto P ER (the positive effect of Bmal to P ER was already present in the original network). Interestingly, this effect was present in the ODE system of [START_REF] Almeida | Transcription-based circadian mechanism controls the duration of molecular clock states in response to signaling inputs[END_REF], as an unbinding term of the PER:CRY complex. As already said, such mass action law kinetic terms are often made implicit by default in the Boolean framework. Here, our method points to the importance of this particular one to avoid unwanted transitions within the attractor. From a modeling point of view, the summary graph thus provides a valuable help in the design of Boolean models of complex, oscillating biological systems.

A tool to assess the general robustness of the attractor

The remaining set of unclassified states U 1 is composed of states where Bmal and P C are off while P ER or CRY can be on. Contrary to U 2 , this set is highly connected to almost every phases of the attractor (see Figure 4). In total, the four edges entering in U 1 amount to 117 transitions. As before, we tried to remove those transitions in order to disconnect U 1 from the rest of the attractor, thus obtaining a final summary graph perfectly matching the order of the phases. However, although showing no more short-cut between phases, this version loses some essential dynamical property. Indeed, the PER and CRY mutants no longer exhibit oscillatory behaviors, suggesting that at least a subset of these transitions are needed for a proper behavior of the model.

Nevertheless, even though those transitions cannot be removed all together, we can still use the summary graph to investigate the model further and analyze each edge entering U 1 separately, in the context of the circadian cycle.

1. The edge V 4 → U 1 is associated with 63 transitions in which the variable P C switches off while P ER and CRY are not both deactivated. From a biological point of view, it means that the PER:CRY complex disappears before PER and CRY have disappeared, indicating a short half-life of the complex. These transitions are to be compared with the V 4 → V 5 transitions where P C switches off after P ER and CRY , suggesting a longer half-life.

In our model the two types of transitions coexist since the relative times of disappearance of the three components, PER, CRY and PER:CRY are not taken explicitly into account. Interestingly, note that the removal of the 63 transitions does not affect the global dynamical properties of the attractor, as the mutants are not affected. 2. The edge V 5 → U 1 corresponds to 9 transitions in the attractor, in which the variable P ER switches on while Bmal is not yet active. Since PER transcription is activated through the CLOCK:BMAL1 complex, these transitions seem rather unrealistic. To investigate further, we removed the transitions and applied the same technique as before, to obtain the modified P ER rule:

P ER = (Bmal ∧ CRY 2) ∨(Bmal ∧ P C) ∨ (DBP ∧ E4)∧(PER ∨ CRY ∨ P C) .
Dynamically, the removal does not affect the main properties of the attractor (the summary graphs of wild type and mutant conditions are similar) however, the inferred logical rule exhibits new and undocumented interactions, namely an auto-activation of PER and a positive effect of CRY on PER. This points to a specific part of the model that will need a closer look in the future. 3. Finally, the edges V 1 → U 1 and V 2 → U 1 correspond to a total of 45 transitions in which the variable Bmal switches to 0 in the absence of P C. With respect to circadian clock events, these transitions describe an early deactivation of Bmal, leading the cycle to bypass the important steps of PER:CRY formation and translocation (phase V 3 ). However, it is the removal of these transitions that directly alter the behaviors of the mutants, suggesting they are necessary to ensure the robustness of the attractor. Their removal leads to a new rule for Bmal: Bmal = (ROR ∧ REV ∧ P C) ∨ (Bmal ∧ (PER ∨ CRY) ∧ PC , that highlights, besides an auto-activation term, a direct positive effect of P ER and CRY on Bmal. Those new interactions generate two positive feedback loops involving Bmal, P ER and CRY that directly interfere with the negative feedback loop between Bmal and P C:

Bmal

) ) P C

P ER, CRY

] ] F F This alteration of the negative loop is directly responsible for the loss of robustness of the model.

The objective of this section was to analyse the "short-cut" transitions in the circadian attractor (see Fig. 4 left). As shown by the three examples above, the summary graph combined with Boolean inference indicates that some of these unwarranted transitions should not be removed and, in fact, appear to play a significant role during the progression of the cycle, namely in the system's response to gene knock-outs.

Conclusions and perspectives

In this article we used the summary graph, a novel tool introduced in [START_REF] Diop | Summarizing complex asynchronous Boolean attractors, application to the analysis of a mammalian cell cycle model[END_REF], to analyze a Boolean model of the mammalian circadian cycle regulation. This tool is well adapted to study biological oscillations, as it provides a rational way to compare a complex Boolean attractor (with hundreds of states) with an oscillating phenomenon. Combined with Boolean inference techniques, it becomes particularly useful in model design as it allows to relate a local dynamical property of the attractor, such as the transition from one phase to the next, to the part of the network's topology directly responsible for it.

After verifying the model's attractor correctly captures the correct succession of phases, the summary graph was further analyzed to provide more insight on the model. In a first step, our analysis indicates that the attractor contains many "short-cuts", that is transitions between states (essentially towards sets U 1 and U 2 ) which may lead to distorted cycles with very short days or very short nights. In a second step, the analysis shows that some of these short-cuts (those passing through set U 2 ) can be removed by a refinement of the logical rules. Namely, the rule of PER should be modified to include the effect of PC dissociation which, for simplicity, was not taken into account in the initial Boolean model (compare Table 2 and equation ( 3)). Finally, in a third step, our analysis suggests that some short-cuts, specifically V 1 → U 1 and V 2 → U 1 , are necessary for a correct performance. These transitions are characterized by an early BMAL1 deactivation, and they are responsible also for generating the short cycle mutants observed in the PER and CRY knock-outs. The transitions through U 1 may thus be necessary to generate circadian cycle robustness in response to perturbations in gene expression.

More generally, the summary graph provides an efficient way to tackle complex qualitative attractors, by testing the effect of specific perturbations on the dynamics. In future works we plan to further investigate the role of the states U 1 in maintaining circadian oscillations, by studying the links between state transitions and the topology of the circadian network.

Fig. 1 .

 1 Fig. 1. Synchronous attractor of the Boolean model (variables are ordered as in Table2). Bold digits indicate the succession of states corresponding to the Bmal/P C phase opposition: in the left column, Bmal is expressed while P C is turned off, the former leading to expression of P ER and CRY (top left); next P C becomes expressed and Bmal is turned off (right column), and CRY and P ER eventually turn off; at the bottom both P C and Bmal are off, before a new cycle begins.

1 |

 1 = 77, |V 2 | = 31, |V 3 | = 31, |V 4 | = 136 and |V 5 | = 30 and there are |U | = 137 unclassified states.

Fig. 3 .

 3 Fig. 3. Summary graph of the attractor on the partition {U, Vi, i = 1 . . . 5} defined in (1), in wild type condition as well as for mutants P ER = 0, CRY = 0 and REV = 0.

Fig. 4 .

 4 Fig. 4. Refined summary graphs. Left: initial attractor's summary graph with U decomposed in (U1, U2); center: summary graph once the six transitions V3 → U2 have been removed; right: summary graph of the new amended model's attractor.

Table 2 .

 2 Logical rules of the Boolean clock model, with 11 variables.

In[START_REF] Tournier | Uncovering operational interactions in genetic networks using asynchronous Boolean dynamics[END_REF], this step is performed by the algorithm REVEAL[START_REF] Liang | Reveal, a general reverse engineering algorithm for inference of genetic network architectures[END_REF], in[START_REF] Tournier | Optimal resource allocation enables mathematical exploration of microbial metabolic configurations[END_REF] it is performed by the algorithm presented in[START_REF] Murakami | Efficient algorithms for dualizing large-scale hypergraphs[END_REF]. Overall, all the inferences performed in the present paper are fast, due to the small dimension (11) of the network. A general discussion about the complexity of the inference problem can be found for instance in[START_REF] Tournier | Optimal resource allocation enables mathematical exploration of microbial metabolic configurations[END_REF].

Acknowledgements

The authors would like to thank Franck Delaunay for many useful discussions on the circadian clock events, as well as for pointers to relevant references.

Supported by the ANR (French agency for research) through project ICycle ANR-16-CE33-0016-01.