M Martens
email: mmartens@fsu.edu

G Franco

N S Dalal

S Bertaina

I Chiorescu

Spin-Orbit Coupling Fluctuations as a Mechanism of Spin Decoherence

We discuss a general framework to address spin decoherence resulting from fluctuations in a spin Hamiltonian. We performed a systematic study on spin decoherence in the compound K6[V15As6O42(D2O)] • 8D2O, using high-field Electron Spin Resonance (ESR). By analyzing the anisotropy of resonance linewidths as a function of orientation, temperature and field, we find that the spin-orbit term is a major decoherence source. The demonstrated mechanism can alter the lifetime of any spin qubit and we discuss how to mitigate it by sample design and field orientation.

INTRODUCTION

In solid-state systems, interactions between electronic spins and their environment are the limiting factor of spin phase lifetime, or decoherence time. Important advances have been recently realized in demonstrating long-lived spin coherence via spin dilution [1][2][3][4][5][6] and isolating a spin in non-magnetic cages [7], for instance. The presence of a lattice can be felt by spins through orbital symmetries and spin-orbit coupling. An isolated free electron has a spin angular momentum associated with a g-factor g e = 2.00232 but in general, spin-orbit coupling changes the g-factor by the admixture of excited orbital states [8] into the ground state. In this Letter, we demonstrate that fluctuations in the spin-orbit interaction can be a significant source of spin decoherence. We present a general theoretical framework to obtain noise spectrum. The method is applied to fluctuations of the long-range dipolar interactions and we observe how the spin-orbit term is modulating the induced decoherence. The model describes spin dilution and thermal excitations effects as well. Experimentally, we analyze shape and orientation anisotropy of ESR linewidths of the molecular compound K 6 [V IV 15 As III 6 O 42 (D 2 O)] • 8D 2 O or V 15 . This system has shown spin coherence at low temperatures [5,9] and interesting out of equilibrium spin dynamics due to phonon bottlenecking [10,11]. However, the details of the spin decoherence are still not fully understood. In the case of diluted or molecular spins, little evidence has been brought up to now on the role of spin-orbit coupling on spin coherence time. This study elucidates this decoherence mechanism and how to mitigate its effect.

FLUCTUATIONS IN SPIN HAMILTONIAN

The V 15 cluster anions form a lattice with trigonal symmetry containing two clusters per unit cell [12]. Individual molecules have fifteen V IV s = 1/2 ions arranged into three layers, two non-planar hexagons sandwiching a triangle (see Fig. 1a). Exchange couplings between the spins in the triangle and hexagons exceed 100 K [13,14] and at low temperatures this spin system can be modeled as a triangle of spins 1/2. The spin Hamiltonian is, as discussed in Supplemental Material[15] (SM) Section I:

H st = H 0 + H J + H DM (1)
where H 0 describes the Zeeman splitting in an external field B 0 , H J is the symmetric exchange term, and H DM is the anti-symmetric Dzyaloshinsky-Moriya (DM) term (see [16] for a detailed formulation). H st eigenvalues are shown in Fig. 1(b) and are used to calculate resonant field positions B res of the ESR spectra through the method of first moments [16]. As shown in Fig. 1(b), the ground state of the total molecular spin S is S = 3/2 for large enough B 0 . In this case, dipolar interactions between total molecular spins in the crystal are described by:

H d = 3µ 0 8π S 2 µ 2 B p;q =p g p (θ)g q (θ) 1 -3 cos 2 φ pq d 3 pq ( 2
)
where µ 0 is the vacuum permeability, µ B is the Bohr magneton, θ is the angle between B 0 and the z axis (z is ⊥ to triangle plane and is also the symmetry c axis of the molecule), d pq is the distance between two molecules located at sites p and q, g p,q (θ) = g 2 a sin 2 θ + g 2 c cos 2 θ 1/2 , g c,a are the g-tensor components parallel and perpendicular to the z axis, φ pq is the angle between S at site p and d pq . Due to local fluctuations of the g-factor, as discussed below, g p and g q are distinct quantities.

The linewidth of ESR signals can be significantly affected by exchange interactions. In V 15 the intramolecular couplings are large and the exchange narrowing effect [17] collapses the (2I + 1) 15 resonances (I = 7/2 for 51 V) into one and it also acts to average out fluctuations related to H st . This leaves fluctuations in H d as being the major contributor to spin decoherence.

There are three possible sources of fluctuation in Eq.( 14), the first being the geometrical factor 1 -3 cos 2 φ pq d -3 pq = R pq (t) since both d pq and φ pq can fluctuate (here, t represents time). This case is described by Bloembergen et al. [18] (Nuclear Magnetic Resonance case) and Kubo and Tomita [17] (ESR case). If R pq (t) fluctuates randomly, its correlation function decays exponentially R(t)R(0) = R 2 + r 2 exp(-t/τ dip ) with a Fourier spectrum:

J R (ν) = 2 π r 2 τ dip 1 + 4π 2 ν 2 τ 2 dip (3)
where R is an average value of the geometric term p =q R pq , r is an average size of R(t)'s fluctuations and the correlation time τ dip is a characteristic of the random motion. This result is described generally by Atherton [START_REF] Atherton | Electron Spin Resonance: theory and applications[END_REF] and can be applied to any stationary random function that is independent of the time origin. The inverse square of the decoherence time T 2 is proportional to J R (ν)dν [17,18]. Therefore, the decoherence rate depends directly on r: 1/T 2 ∝ r.

Another fluctuation source comes from thermal excitations to different S Z states of S, where Z axis is B 0 , which defines the second moment of a resonance line [START_REF] Pryce | [END_REF][21][22] (potential fluctuations between different spin states in low fields has been studied as well [9]):

S Z (t) 2 S Z (0) 2 = S 4
Z + KU (T ) exp(-t/τ s ), where τ s is the thermal correlation time and KU (T ) a term studied by Kambe and Usui [21]. It is shown that the fluctuations Fourier spectrum is proportional to a temperature dependent factor:

KU (T ) =< S 2 Z > T -< S Z > 2 T = S 2 d dy B s (y) (4)
where B s (y) is the Brillouin function, y = T Z S/T , T Z = hf 0 /k B (f 0 is the microwave excitation frequency), and S = 3/2 is the total spin state. KU (T ) has thus a similar role to r 2 in Eq. ( 3). This formulation is valid above the ordering temperature which is ∼0.01 K [23] for V 15 .

The dipolar term H d serves as an excellent platform to study fluctuations of g(θ). Its value away from g e is due to the spin-orbit interaction and it is given by [8]:

g = g e I -2λΛ (5)
where g is the g-tensor (diagonal [g a , g a , g c ] for V 15 ), I is the unit matrix, λ is the spin-orbit coupling constant and Λ is a tensor defined in terms of the matrix elements of the orbital angular momentum L. In general terms, Λ is the coupling between the ground and excited orbitals divided by their energy separation. Relative fluctuations with an average size ξ = δ(λΛ)/(λΛ) (assumed isotropic) can be induced by crystal and molecular vibrations. In particular, Raman measurements on V 15 [24] discussed below, show a broad distribution of the vibration modes. Fluctuations of excited orbitals and thus of Λ can generate broad virtual transitions since those orbitals are mixed with the ground orbital state. The resultant fluctuation in the g-factor can be written as: Assuming g(θ) is a stationary function with small temporal random fluctuations and that magnetic and orbital fluctuations are uncorrelated in first approximation, the correlation function of a fluctuating H d (t) is:

δg(θ) = ξ (g(θ) -g e ) . (6)
G d (t) = H d (t)H d (0) = α 2 g(t)g(0) 2 S 2 Z (t)S 2 Z (0) R(t)R(0) (7)
where g(t)g(0) = g(θ) 2 + (δg(θ)) 2 exp(-t/τ g ), τ g is the correlation time of g-factor fluctuations, and α = 3µ0µ 2 B 8π . A corresponding J d (ν) gives the Fourier spectrum of the fluctuations, as in Eq. 3. G d (t) can be written as the sum of four terms (see SM Section II for details): G 0 which is a constant , G g (t) ∝ g(θ) 4 , G δ (t) which is temperature indepedent and G T (t) which is temperature dependent.

In absence of g-factor fluctuations, the resulting Fourier spectrum is defined only by G g (t) and for neglible r (less important in solids at low temperatures) the term

G g (t) ≈ α 2 g(θ) 4 R 2 KU (T )e -t
τs is as in [21]. A temperature dependence of the linewidth ∝ KU (T ) is similar to observations done with Fe 8 [25][26][27], nitrogen-vacancy color centers in diamond [28] while other studies seem to confirm the proportionality to the g-factor [29,30]. If the g-value does fluctuate then all three terms G g,δ,T represent sources of decoherence, with G δ + G T given by:

G δ (t) + G T (t) ≈ α 2 R 2 S 4 + KU (T ) × × 2g(θ) 2 δg 2 (θ)e -t τg + δg 4 (θ)e -2t τg . (8)
Because 1/T 2 2 is ∝ J d (ν)dν, an important consequence is that one can combine different decoherence sources by summing their effect (each term i) as follows:

1 T 2 2 ≈ i 1 T 2 2i
, similar to the well-known fact that the sum of uncorrelated variances is equal to the total variance. Additionally, the weight of each term in the sum depends on, or can be tuned with, the field angle θ through g and δg. Here we show that for V 15 , the anisotropy of the decoherence time is explained by fluctuations δg, as shown in Eq. 8, amplified by spin thermal fluctuations KU (T ).

EXPERIMENTAL DATA

Continuous-wave ESR measurements at 120, 241, and 336 GHz are performed using the quasioptical superheterodyne spectrometer at the National High Magnetic Field Laboratory [31,32], with a sweepable 12.5 T superconducting magnet (homogeneity of 10 -5 over 1 cm 3 ). Sample temperature can be varied from room temperature down to 2.5 K. A single crystal of regular shape (as in [16]) of volume < ∼ 0.1 mm 3 was positioned on a rotating stage allowing for continuous change of the angle θ between B 0 and the c axis of the molecule following the procedure described in [16]. The homogeneity of the magnet compared to the size of the crystal allows us to ignore B 0 as a source of broadening. The applied fields are above 4 T, past the crossing of the S = 1/2 doublet and S = 3/2 quartet, such that the ground state of the system is in the S = 3/2 quartet (see Fig. 1(b)).

ESR spectra at temperatures T = 4-60 K for B 0 and ⊥ to c-axis (θ = 0 • , 90 • respectively) show a Lorentzian (homogenous) lineshape. Representative spectra with Lorentzian and Gaussian fits are shown in Fig. 2(a) for comparison. The temperature dependence of the linewidth is shown in Fig. 2(c) for three microwave frequencies f 0 . Compared to measurements made at lower fields [33], where the ground state is in the S = 1/2 doublet, the linewidths are ∼ 10 times narrower. Plotted is the full width at half maximum (FWHM) of the Lorentzian fits vs T /f 0 to underline that the temperature dependent mechanism of decoherence in the system qualitatively follows the temperature behavior predicted by KU (T ) plus a temperature independent contribution, essential in the low T limit [Eq. ( 8)]. Note the absence of linewidth increase with frequency (or field) which excludes a static distribution of the g-factor. Additionally, there is no hyperfine structure visible in the spectra (exchange narrowing) since the 3d electrons of V interact with the nuclei of several other V ions due to the large exchange couplings (∼ 10 2 K) within the molecule. Due to these properties of the measured line width and shape, we can estimate T 2 to be the inverse of the FWHM.

There are two distinct curves in Fig. 2(c), dependent on the orientation of B 0 . To probe this orientation dependence, the linewidth is measured as a function of θ, see Fig. ( B 0 c axis) while the largest occurs at θ = 90 • ( B 0 in the triangle plane). This implies more decoherence the further g(θ) gets from g e since g(0 • ) = g c ≈ 1.98 and g(90 • ) = g a ≈ 1.95. The fact that the width is largest (smallest) when g(θ) is minimum (maximum) rules out exchange narrowing being the cause of this anisotropy since it would require a linewidth ∝ (1 + cos 2 θ) [17], as in the case of CsCuCl 3 [34]. To rule out an angular effect of the dipolar field distribution, we measured the FWHM(θ) on a sample of irregular shape at 240 GHz and 60 K. Although the shape-dependent coefficients R and r must be different, the same behavior is observed as in Fig. 2b (see SM Section III B for details). We can thus focus on the terms G g,δ,T as source of fluctuations.

However, the term G g (t) is ∝ g(θ) 4 , in clear contrast with the observation that the FWHM and g(θ) have opposite angular dependences (see Fig. 2b and also SM Section IIIA). This is the essential property of the V 15 compound, which makes it particularly suitable to study the effect of spin-orbit fluctuations. Therefore, this opposite angular behavior provides evidence that δg(θ) = 0 and the terms G δ,T (t) must be considered while G g (t) can be discarded. One can argue that geometrical fluctuations in solids at low temperatures are very small (r ≪ R) and lattice fluctuations are mostly influencing the relaxation time T 1 (τ s ≫ τ g ) making G g ≈ constant at the time scale of the decoherence time.

Since [17,18] the linewidth square can be modeled by the following fit function (see SM Section II for details):

G d (t) ≈ G δ (t) + G T (t) and 1/T 2 2 ∝ J d (ν)
∆ 2 = S 4 + KU (T ) × × 2a 2 g(θ) 2 (g(θ) -g e ) 2 + A 2 (g(θ) -g e ) 4 (9)
where A and a are fit parameters. The procedure is detailed in SM Section IV; it allows to calculate the angular dependence FWHM(θ) by using only two data points, ∆(0 • ) and ∆(90 • ), as shown in Fig. 2

(b) (dashed lines).

To analyze the temperature dependence of the linewidth shown in Fig. 2(c), we solve for A and a at all available temperatures and frequencies (see SM Section IV B for details). Above 10-20 K, the values stabilize at A ∼100 GHz and a ∼3.2 GHz. At lower temperatures, the values decrease by almost half, indicating a small decrease in ξ and/or a slowing down in the fluctuations time τ g . These temperature trends A(T ) and a(T ) are estimated by an exponential saturation (see SM Fig. 4 for details), with decay constants of 3.6 K and 11 K for A(T ) and a(T ) respectively. With no other adjustments, the calculated linewidth is in very good agreement with the experimental data, as shown with dashed line in Fig. 2(c). On the low end of T /f 0 one observe a residual value of the linewidth, which includes the effects of other decoherence sources (such as the nuclear spin bath [35]), although it can be well described by Eq. 9.

The outcome of the fit procedure can be used to estimate the size of spin-orbit fluctuations (see SM Section IV C) leading to an order of magnitude for ξ ∼ 10 -2 . This corresponds to a fluctuation δg/g ∼ 10 -4 , too small to result in directly measurable fluctuations of the Zeeman splitting. Note that for V 15 , a large spin-orbit fluctuation is supported by previous Raman measurements [24] showing a very broad signal in the region of ∼ 500 cm -1 corresponding to vibrations of oxygen bridges between V ions. The observed broad distribution of the modes can induce very fast virtual transitions to excited coupled states and, as a consequence, spin decoherence.

CONCLUSION

Our study provides insight on how to mitigate the effects of spin-orbit fluctuations. It is evident from Eqs. ( 6) and ( 7) that the g-tensor should be as close as possible to g e . In molecular compounds this can be achieved by engineering the ligands type since local symmetry affects the diagonal values of the g-tensor of a magnetic ion. Aside from material design by chemical methods, J d (ν) can be minimized by applying the magnetic field at a specific angle θ. For V 15 , this would be θ = 0 for which the decoherence time reaches several nanoseconds. This time can reach ∼ 400 ns by reducing R in J d (ν) via dilution in liquid state, thus allowing the observation of Rabi oscillations and spin-echoes [5]. The methodology presented here can be important for the diluted spin systems as well, since long range interactions are still present and can carry modulations due to g-factor fluctuations. Potential examples are transition metals such as Cr 5+ :K 3 NbO 8 or some lanthanide monomers doped into insulating lattice such as Hf 3+ :LuPO 4 [36] or La 2+ :CaF 2 [37]. The results extend to any solid-state system where spin-orbit coupling leads to quantum effects, independent of system dimensionality.

SUPPLEMENTAL MATERIAL: SPIN-ORBIT COUPLING FLUCTUATIONS AS A MECHANISM OF SPIN DECOHERENCE SINGLE MOLECULE SPIN HAMILTONIAN

The molecular compound K 6 [V IV 15 As III 6 O 42 (D 2 O)] • 8D 2 O or V 15 contains fifteen V +4 spins, each with a spins 1/2, shown in blue in Figure 1(a) of the article. The spins are coupled by large antiferromagnetic (AF) couplings in exces of 100 K [13,14]. At low temperatures, the molecule can be modeled as an effective frustrated triangle of spins 1/2. The effective antiferromagnetic coupling is ≈ 2.51 K as estimated by magnetic measurements at low temperature [23]. The generalized three spin Hamiltonian H st is given in the article as:

H st = H 0 + H J + H DM . ( 10
)
where H 0 describes the Zeeman splitting in an external field B 0 , H J is the symmetric exchange term, and H DM is the anti-symmetric Dzyaloshinsky-Moriya (DM) term.

In the study presented here, the focus lies on the fluctuations of the dipolar term. However, for completness of the study, a summary of the properties of the static term H st is presented in this section following our previous study, Ref. [16].

The Zeeman term is given by:

H 0 = g a µ B B △ S △ + g c µ B B z S z (11)
where g a,c are the values of the g-factor in the triangle's plane and along the molecule's c axis (z axis in Fig. 1a) respectively; B △,z are the planar and vertical vector components of the external magnetic field B 0 , with B △ = B 0 sin θ and B z = B 0 cos θ, respectively; θ is the angle between B 0 and the z-axis; S △ = S x + S y = i=1,2,3 ( S (i)

x + S (i) y ), and S z = i=1,2,3 S

z . The exchange term is given by:

H J = <i,j> [J (ij) z S (i) z S (j) z +J (ij) t (S (i) x S (j)
x +S (i) y S (j) y )] ( 12)

and the DM term by:

H DM = <i,j> D (ij) ( S (i) × S (j) ) ( 13
)
where one sums over the (i, j) pairs (1, 2), (2, 3), (3, 1) of the triangle spins S (1),( 2), (3) . The H J term describes the effective antiferromagnetic coupling between the triangle's corners which, in principle, can be anisotropic (J z = J t ) and differ between spin-pairs in the triangle (J

(12) z = J (23) z
for example). In the H DM term, the individual vectors

D (ij) = D (ij) x + D (ij) y + D (ij) z
correspond to the interaction between S (i) and S (j) , and are related to the local frames (a, b, and c) shown in Figure 3(a). It = J 0 = -2.45 K for all (ij), g a = 1.95 and g c = 1.98. In this simple case, the ground doublets are degenerate and are separated from the quartet at zero field by 3J 0 /2. The crossover of the ground state from the doublet to the quartet is known to happen at 2.8 T from measurements of the magnetization of a single V 15 crystal, [23] placing a constraint on the values of J

(ij) z and J (ij) t .
The parameters of H st can be identified by means of high-field spin resonance spectroscopy as described in a previous study [16]. The excitation frequency of 120 GHz corresponds to fields upwards of 4 T, where the quartet represents the ground state. The allowed transitions are shown in Figure 3(b) with dashed vertical lines, where the blue line between the crosses represents the four allowed transitions within the degenerate doublets. The field positions of these transitions depend on the values of all H st terms and each transition has a certain effect on the overall field location of the resonance signal. By using numerical diagonalization and the first moment method to average the transitions, information on the spin Hamiltonian parameters can be extracted. Such approach is very general and can be applied to solve a non-linear relationship between magnetic field and resonance frequency. The conclusion of the study can be summarized as follows. Any isotropy in the AF exchange coupling J or generated by the DM term is not sufficiently large to create a spectroscopic signature. The data together with the numerical calculations of the transition probabilities indicate that the doublet transitions, shown in blue arrows in Figure 3(b), don't have any noticeable influence on high-field spin resonance experiments.

SPECTRAL DENSITY OF DIPOLAR FLUCTUATIONS

The interaction between V 15 neighboring molecules is of dipolar nature and described by Eq. ( 2) of the main article:

H d = 3µ 0 8π S 2 µ 2 B p;q =p g p (θ)g q (θ) 1 -3 cos 2 φ pq d 3 pq ( 14
)
with µ 0 the vacuum permeability, µ B the Bohr magneton, d pq the distance between two molecules located at sites p and q, g p,q (θ) = g 2 a sin 2 θ + g 2 c cos 2 θ 1/2 , φ pq the angle between S at site p and d pq . Random fluctuations in H d generate the observed linewidth and thus limit the spin coherence lifetime.

To find the spectral density of dipolar fluctuations one has to analyse the time dependence of the correlation function, defined in the main article as

G d (t) = H d (t)H d (0) = α 2 g(t)g(0) 2 S 2 (t)S 2 (0) R(t)R(0) : G d (t) = α 2 g(θ) 2 + δg 2 (θ)e -t τg 2 × S 4 + KU (T )e -t τs R 2 + r 2 e -t τ dip = α 2 g(θ) 4 + 2g(θ) 2 δg 2 (θ)e -t τg + δg 4 (θ)e -2t τg × S 4 + KU (T )e -t τs R 2 + r 2 e -t τ dip (15)
where all notations are defined in the main article. Thus, G d (t) is a twelve term sum which can be split as follows:

G 0 = α 2 g(θ) 4 S 4 R 2 , G g (t) = α 2 g(θ) 4 R 2 KU (T )e -t τs + + α 2 g(θ) 4 r 2 e -t τ dip S 4 + KU (T )e -t τs , G δ (t) = α 2 S 4 2g(θ) 2 δg 2 (θ)e -t τg + δg 4 (θ)e -2t τg × R 2 + r 2 e -t τ dip , G T (t) = α 2 KU (T )e -t τs 2g(θ) 2 δg 2 (θ)e -t τg + δg 4 (θ)e -2t τg × R 2 + r 2 e -t τ dip . ( 16
)
As discussed in the main article, one can follow Atherton [START_REF] Atherton | Electron Spin Resonance: theory and applications[END_REF] to obtain the Fourier spectrum corresponding to each term above. The first term, G 0 , is time indepedent and thus will not count. Also, given the experimental data for the particular case of V 15 , the term G g must not be significant. This is a result of the opposite angular dependence between measured g(θ) and signal linewidth, an analysis given in Section . This condition imposes that r 2 is very small, meaning that crystal lattice fluctuations at low temperatures are not of significance for the decoherence rate. Such process is significant for the relaxation time T 1 , related to τ s and usually larger than the decoherence time. Consequently, τ s ≫ τ g , τ dip and therefore the first term in G g (t) is almost constant, at the scale of the decoherence time. Such condition is valid in a large range of temperatures and fields, leaving only two dominant terms, G δ and G T .

The essential difference between G δ and G T is that the former is temperature independent while the later describes the measured temperature dependence (see Fig. 2b). At the low end of the temperature range, one achieve a minimization of the continuous wave linewidth. This can be attributed to the temperature independent term G δ plus any other decoherence sources, such as the nuclear spin bath [35].

Consequently, the correlation function can be written as:

G d (t) ≈ G δ (t) + G T (t) ≈ α 2 S 4 + KU (T ) R 2 × 2g(θ) 2 δg 2 (θ)e -t τg + δg 4 (θ)e -2t τg , (17)
and the resulting Fourier spectrum is:

J d (ν) = 2 π α 2 S 4 + KU (T ) R 2 × 2g(θ) 2 δg 2 (θ) τ g 1 + 4π 2 τ 2 g ν 2 + δg 4 (θ) τ g /2 1 + π 2 τ 2 g ν 2 . ( 18
)
Since the integral of J d (ν) represents the square of the decoherence rate, 1/T 2 2 = J d (ν)dν, one can conclude that the appropiate fit function for the spectroscopy linewidth is given by:

∆ 2 = 2 π α 2 S 4 + KU (T ) R 2 × 2g(θ) 2 ξ 2 (g(θ) -g e ) 2 τ g 1 + 4π 2 τ 2 g ν 2 dν + + 2 π α 2 S 4 + KU (T ) R 2 × ξ 4 (g(θ) -g e ) 4 τ g /2 1 + π 2 τ 2 g ν 2 dν , (19)
where ξ = δg(θ)/(g(θ)-g e ) and g e = 2.00232. This leads to a fit function of the experimental linewidth ∆:

∆ 2 = S 4 + KU (T ) × 2a 2 g(θ) 2 (g(θ) -g e ) 2 + A 2 (g(θ) -g e ) 4 , (20)
where a and A are fit parameters.

ANGULAR DEPENDENCE OF THE LINEWIDTH/SPECTRAL DENSITY

The g 4 term is negligible

When discussing the correlation function G d (t) above, an approximation was made, namely that the G g (t) term must not be dominating because it contains the factor g(θ) 4 . The V 15 molecular compound has a very particular propriety, namely the behavior of the linewidth when the angle θ is changed continously. Our measurements (see Fig. 2b of the main article) show an opposite angular dependence between linewidth and g(θ). At the extreme points θ = 0 • and 90 • a small/large linewidth is obtained when the g-factor is the largest/smallest, with g(0 • ) = g c ≈ 1.98 and g(90 • ) = g a ≈ 1.95. This is somewhat counter-intuitive since a large g-factor means a larger magnetic moment and thus potentially larger dipolar fluctuations. This unusual behavior can not be described by the term g 4 (θ) which corresponds to a traditional interpretation of dipolar fluctuations: the larger the value of the dipoles (∝ g, S) the larger the value of H d fluctuations. Instead, the linewidth is in sync with the angular dependence of δ(θ), that is the size of the shift of the g-factor away from the free electron value g e . This shift of g is caused solely by spin-orbit coupling [8]: g = g e I -2λΛ where g is the g-tensor (diagonal [g a , g a , g c ] for V 15 ), I is the unit matrix, λ is the spin-orbit coupling constant and Λ is a tensor defined in terms of the matrix elements of the orbital angular momentum L:

Λ ij = n =0 0|L i |n n|L j |0 (E n -E 0 ) (21)
with i, j = x, y, z. The spin-orbit coupling to excited orbitals is defining the value of the g-factor and vibrational modes will induce fluctuations in the term λΛ. As a consequence to the opposite angular behavior of the linewidth and the g-factor, we conclude that δg fluctuations with a characteristic time τ g modulate the dipolar fluctuations. This makes the terms in δg 2 and δg 4 , rather than g 4 , to dominate the spectral density J d (ν) shown above.

Dipolar Hamiltonian: the geometric factor

As a consequence of the long range dipolar interactions, each molecule will see a local field summing the individual dipolar fields generated by the other molecules. Therefore, the local field has a certain total value and also a standard deviation due to the spread in values of the individual fields. By means of numerical simulation, one can analyse the dependence on θ of the spread in dipolar fields and compare it to the observed angular dependence of the experimental resonance linewidths. As a reminder, the main observation in the case of V 15 is that the linewidth is small/large when the g-factor is large (θ = 0 • ) / small (θ = 90 • ).

The angular dependence of the distribution of fields is calculated in the middle of an ensemble having 22 molecules on each crystallographic axis, that is a total number of 22 × 22 × 22 = 10648 sites. The size is chosen such that dipolar couplings due to molecules outside of this interaction volume can be neglected. The ensemble is shown in the insert of Figure 4. The applied field B 0 is large enough to polarize the spins. The total field along B 0 is given by B dip = q b dip q where b dip q is the individual field generated by a molecule at site q in the middle of the crystal. The result is plotted with a continous red line, while the spread σ 2 dip = q (b dip q ) 2 -( q b dip q ) 2 is shown with a dashed blue line as a function of θ. One note that σ dip is significantly smaller than the large linewidths observed experimentally (about an order of magnitude smaller) and also that the simulated angular dependence is not correlated with the measured one. The simulated B dip shows the typical behavior ∝ 1 -3 cos 2 θ and both B dip and σ dip do not show the experimentally observed behavior of the linewidth (minimum/maximum values at θ = 0 • /90 • respectively).

This shows that fluctuations in the geometrical part of H d can not explain the observed angular dependence of the linewidth.

Another way to verify the above statement is to measure a sample of an irregular shape as shown in the insert of Figure 5. Given that dipolar interactions are of long range nature, the shape-dependent coefficients R and r could completely change the angular dependence of the linewidth. The data in such case can be analyzed and compared with the case of a well defined crystal, with a shape similar to the insert of Figure 4 and as the one studied in the main article. The linewidth measured at T = 60 K and excitation frequency 240 GHz is shown with red dots in Figure 5. Despite of an irregular shape of the crystal, the main property of the linewidth is conserved, namely that its angular dependence is opposite to that of the g-factor. We are thus compelled to relate the linewidth to fluctuations of terms in H d containing the g-factor .

FITTING PROCEDURE Field-frequency conversion

The measured field linewidths can be converted into frequency units through exact diagonalization of the three-spin Hamiltonian H st , since the linewidth corresponds to a broadening of its eigenvalues. This is needed since the resonance frequency and field are not in a linear relantionship due to zero field splittings. The method presented below is general and can be applied to any other molecular system. We calculate the minimum and maximum excitation frequencies f min,max leading to resonance fields B max,min using the first moment method [16]. Their difference is:

∆ = f max -f min . (22)
By measuring ∆(0 • ) and ∆(90 • ), one can construct a system of two equations with two unknowns (see Eq. 20) to solve for A and a. This allows to calculate the curve ∆(θ) which is converted back to field units by solving the static Hamiltonian: for two frequencies f ′ min,max (θ) = f 0 ± ∆(θ)/2, we calculate the corresponding resonance fields B ′ max,min (θ). The difference B ′ max -B ′ min is the fitted Full Width at Half Max (FWHM) as a function of θ. Thus, by using only two data points (linewidths at θ = 0 • and 90 • ) one can simulate the full angular dependence of the linewidth, for any non-linear relantionship between applied field and resonance frequency.

Frequency and temperature studies

The protocol described in the previous section can be applied to find the fit parameters A and a in the temperature interval 4-60 K for the three frequencies used in our study: 120 GHz, 241 GHz and 336 GHz. The obtained values for the fit parameters A and a are shown with scattered dots in Figure 6, in panels (a) and (b) respectively: red squares for 120 GHz, black triangles for 241 GHz and blue circles for 336 GHz. The temperature dependence of their averages A v (T ) and a v (T ) is shown with a smoothed dashed green line lying inside an uncertainty area bounded by two purple dotted, smoothed lines A v ± σ A and a v ± σ a where σ A,a are the standard deviation of the means.

For both fit parameters one can note an increase with the temperature up ∼ 10 -20 K when they reach more stable values. Since A and a reflect the magnitude of the characteristic time τ g and ξ, one can conclude that the spin-orbit fluctuations are slightly less effective at low temperatures. It is important to take this aspect into account when using our model to fit the T/F master curve shown in Fig. 2(c) of the main article. In Figure 6, the black continous lines A(T ) and a(T ) show exponential approximations of the temperature trends, with characteristic decay constants of 3.6 K and 11 K for A and a, respectively. Different decay constants can be linked to different power dependence on ξ for A and a, as shown by Eq. 23. These exponential decays are roughly within the limits of A and a uncertainty bands and are optimized to properly fit the master curves in Fig. 2(c) of the main article (shown with dashed lines).

Discussion of model variables

The values of parameters A and a can be used to further analyze the variables used by the model given in Eq. 19. The integrals over the Lorentzian spectra are from 0 to ∆ (see Ref. [18]) which leads to:

a 2 =
2/π 2π α 2 ξ 2 R 2 arctan(2πτ g ∆)

A 2 = 2/π π α 2 ξ 4 R 2 arctan(πτ g ∆). ( 23
)
The g 2 δg 2 term, which defines the value of a, is the leading term of the fit procedure and fully describes the angular behavior of the linewidth; the δg 4 term provides only a marginal improvement of the fit. Consequently, the parameter A is particularly susceptible to a large uncertainty due to the small value of the δg 4 term: in order to have a meaningful contribution to the fitting procedure, A has to be quite large. Therefore, taking the ratio A/a as an estimation of ξ can lead to unrealistic values for the size of the g-factor fluctuations.

Instead, one can use a to asses the value of ξ, after evaluating R. As a reminder, the coefficient R is an average of p =q R pq (as defined in the main article) and it depends on the spin density and sample shape. An estimation of R can be made by fixing one site, as in the simulations presented in Section III B, that is R = N q R pq with p fixed in the center of the crystal and N = 10, 648 is the total number of sites in the interaction volume. The dipolar field B dip , shown in Fig. 4 with a red continous line, can be estimated as B dip ∼ µ0 4π gµ B SR/N , which for g = 2, S = 3/2 and B dip ≈ 0.5 mT, leads to R ∼ 2000 nm -3 .

Since the linewidth of the homogenous ESR signal in V 15 is quite large, compared to other spin systems proposed as qubits, one can assume that τ g ∆ > 1, which is sufficiently large to approximate the arctan function with π/2 and thus a ∼ αξR. Assuming a value of a of ∼ 2 GHz, one get an estimation of the order of magnitude to be ξR ∼ 100 nm -3 and ξ ∼ 0.01 -0.1 as an order of magnitude.

FIG. 1 .

 1 FIG. 1. (Color online) (a) Ball-and-stick representation of V15 (V ions in blue). The x axis is along one side of the triangle while the z axis is perpendicular to the triangle plane and represents the c axis of the crystal unit cell. (b) Level diagram of the three spin model in field ||z, with positions of the three experimental frequencies shown. Dashed lines show the S = 1/2 doublets with the red dashed arrows indicating those transitions. Lines show the S = 3/2 quartet with blue arrows indicating the transitions; the resonance fields are averaged in the first moment calculation of Bres at a given frequency.

FIG. 2 .

 2 FIG. 2. (Color online) (a)Typical measurements of the derivative of the absorption χ ′′ at 120 GHz, 241 GHz and 336 GHz with derivative of Gaussian (blue dotted line) and Lorentzian fits (red dashed line). (b) FWHM of Lorentzian fits as a function of field angle θ measured at three frequencies: 336 GHz (blue triangles), 241 GHz (red circles) and 120 GHz (black squares) The dashed lines are calculated widths as a function of θ; the agreement shows the predicted correlation between decoherence rate and ge -g(θ). In contrast, the green line (right axis) shows the opposite angular behavior of calculated g(θ), leading to Gg ≪ G δ + GT (see text). (c) FWHM of Lorentzian fits vs temperature/frequency for the 3 studied frequencies. Dashed lines are calculated FWHM(T) for θ = 0 • and 90 • .

2 FIG. 3 .

 23 FIG. 3. (color online) (a) Diagram showing the local a, b, c axes defining the DM vectors in the three spin model. The spins are located in the triangle corners. The applied field B0 makes an angle θ with the z-axis ⊥ to the triangle plane. (b) Energy eigenstates of the three spin model with positions of allowed transitions at 120 GHz shown. Red dashed lines between red open circles indicate transitions inside the S = 3/2 quartet while the blue dashed line between the blue crosses indicate the four allowed transitions within the degenerate S = 1/2 doublets. The zero-field splitting between the quartet and doublets is 3J0/2.

=

  D z . Diagonalizing H st yields an energy structure that is made up of an excited quartet S = 3/2 and two ground doublets S = 1/2 where the zero-field splitting between the quartet and doublets is determined by the value of the antiferromagnetic coupling. This is shown in Figure 3(b) for the case where there is no DM interaction, J

FIG. 4 .

 4 FIG. 4. (color online) Dipolar field (continous line) and σ dip (dashed line) as a function of θ calculated in the middle of an ensemble of 10648 molecular spins, shown in the insert. The B dip has a typical ∝ 1 -3 cos 2 θ behavior.

FIG. 5 .

 5 FIG. 5. (color online) Experimental linewidth (red dots) measured at 60 K and 240 GHz for a crystal of irregular shape (see insert) as a function of θ. Despite the irregular shape, the linewidth shows the same overal behavior with a minimum/maximum at θ = 0 • /90 • , respectively.

FIG. 6 .

 6 FIG. 6. (color online) Fit parameters A and a as a function of temperature, shown in panels (a) and (b) respectively, for three excitation frequencies. Their average is shown with a smoothed dashed green line; two purple dotted, smoothed lines are bounding the uncertainty area defined as the average ± the standard deviation of the mean. The black continous line shows an exponential approximation of the temperature trend, used to calculate the temperature / frequency fit (see Fig. 2(c) of the main article).

ACKNOWLEDGEMENTS

We wish to acknowledge David Zipse and Vasanth Ramachandran for their help in growing V 15 crystals. This work was supported by NSF Grant No. DMR-1206267 and CNRS-PICS CoDyLow. The NHMFL is supported by Cooperative Agreement Grant No. DMR-1157490 and the state of Florida. The experiments were performed at the EPR facilities of the NHMFL, and we thank dr. Johan van Tol for his support.


Ovidiu Radulescu

Sandrine Lagarrigue

Anne Siegel

Philippe Veber

Michel Le Borgne

Topology and static response of interaction networks in molecular biology

Keywords: Systems biology, gene and metabolic networks, steady state shifts, Dirichlet to Neumann map

Introduction

Network-based representations are widely used to describe gene regulation or metabolic pathways at cellular level; graph theoretical properties of these biological networks, dynamical properties, as well as their inference from experimental data, were the subject of intensive studies [START_REF] Kaminski | Practical approaches to analyzing results of microarray experiments[END_REF][START_REF] Yamanishi | Protein network inference from multiple genomic data: a supervised approach[END_REF][START_REF] Wagner | How to reconstruct a large genetic network from n gene perturbations in fewer than n 2 easy steps[END_REF][START_REF] Kholodenko | Untangling the wires: a strategy to trace functional interactions in signaling and gene networks[END_REF][START_REF] Thieffry | Qualitative analysis of gene networks: Towards the delineation of trans-regulatory modules[END_REF][START_REF] De | Modeling and simulation of genetic regulatory systems: A literature review[END_REF][START_REF] De Jong | Qualitative simulation of genetic regulatory networks using piecewise-linear models[END_REF]. However, network inference techniques need huge amounts of data which are not always available or often have poor accuracy. Similarly, dynamical studies need accurate data.

In this paper, our goal is different. Instead of using data and gene perturbations for building networks from scratch, we develop mathematical techniques in order to refine the analysis of incomplete models or to compare models and data. Our mathematical results connect network topology and the response to steady state shift experiments. Steady state shift experiments are useful tools in chemistry allowing in principle to recover the reaction mechanisms [START_REF] Chevalier | Toward a systematic determination of complex reaction mechanisms[END_REF]. We argue that similar approaches are well adapted to differential microarray experiments which compare gene expressions between two different states [START_REF] Kaminski | Practical approaches to analyzing results of microarray experiments[END_REF]. † Author for correspondence ovidiu.radulescu@univ-rennes1.fr.

A theory of steady state shifts is in fact a theory of static response, linear when the shifts are small and non-linear when the shifts are big. Theories of linear response were developed in various contexts such as condensed matter and in statistical physics [START_REF] Kubo | Statistical Physics II[END_REF], electrical, physical and mechanical systems [START_REF] Macfarlane | Dynamical Sytems Models[END_REF], chemical and biochemical systems [START_REF] Oster | Network thermodynamics: dynamic modelling of biophysical systems[END_REF][START_REF] Perelson | Chemical reaction dynamics. part ii: Reaction networks[END_REF], complex fluids [START_REF] Larson | Constitutive Equations for Polymer Melts and Solutions[END_REF], metabolic networks [START_REF] Kacser | The control of flux[END_REF][START_REF] Heinrich | A linear steady-state treatment of enzymatic chains. general properties, and effector strength[END_REF] and gene networks [START_REF] Kholodenko | Untangling the wires: a strategy to trace functional interactions in signaling and gene networks[END_REF][START_REF] Vlad | Response experiments for nonlinear systems with application to reaction kinetics and genetics[END_REF]. The basic quantities in such theories are the susceptivities, representing derivatives (generically functional derivatives) of outputs with respect to inputs. Susceptivities (and also response in general) are obtained from constitutive equations that must be compatible with thermodynamics [START_REF] Oster | Network thermodynamics: dynamic modelling of biophysical systems[END_REF][START_REF] Grmela | Complex fluids subjected to external influences[END_REF]. Nonlinear response theories that apply to large perturbations are generally more difficult to handle [START_REF] Larson | Constitutive Equations for Polymer Melts and Solutions[END_REF][START_REF] Grmela | Complex fluids subjected to external influences[END_REF]. [START_REF] Vlad | Response experiments for nonlinear systems with application to reaction kinetics and genetics[END_REF] showed that under special "neutrality" conditions the response of reaction-diffusion systems to special perturbations consisting in changing the abundance of some marked individuals is linear even at large perturbations.

In the static case, structural stability [START_REF] Ruelle | Elements of differentiable dynamics and bifurcation theory[END_REF][START_REF] Smale | Mathematics of Time[END_REF] suggests a simple approach to nonlinear response. Like in classical thermodynamics [START_REF] Callen | Thermodynamics and a Introduction to Thermostatics[END_REF], finite static response can be obtained by integrating the (state dependent) differential forms whose coefficients are the susceptivities provided that the integration path does not contain singularities (phase transitions).

Susceptivities depend on the network topology and therefore can be used in reverse-engineering. This method is complementary to already well developed correlation analysis of gene and metabolic networks [START_REF] Kaminski | Practical approaches to analyzing results of microarray experiments[END_REF][START_REF] Yamanishi | Protein network inference from multiple genomic data: a supervised approach[END_REF] or of time series in chemical reaction systems (A. Arkin and Ross, 1995). [START_REF] Vlad | Response experiments for nonlinear systems with application to reaction kinetics and genetics[END_REF] showed that under the "neutrality" condition the dynamical susceptivity of a system of chemical reactions with respect to abundance perturbations is the exponential of a connectivity matrix, which contains information about the chemical reaction mechanism. It is known [START_REF] Chevalier | Toward a systematic determination of complex reaction mechanisms[END_REF] that small concentration shifts measurements allow to calculate the Jacobian of a chemical kinetics evolution equations. Very related to this is the work of Sontag [START_REF] Kholodenko | Untangling the wires: a strategy to trace functional interactions in signaling and gene networks[END_REF] and the older Kacser-Burns connectivity property from metabolic networks [START_REF] Cornish-Bowden | Fundamentals of enzyme kinetics[END_REF], meaning that the matrix of elasticities (the Jacobian) is the inverse of the matrix of control coefficients (static susceptivities). The dynamical Jacobian is also important in gene networks dynamics. It defines the interaction graph [START_REF] Soulé | Graphic requirements for multistationarity[END_REF], which is the natural mathematical candidate to represent the biologist's gene regulation graph [START_REF] Thomas | On the relation between the logical structure of systems and their ability to generate multiple steady states or sustained oscillations[END_REF]. Systems of chemical reactions allow for other graphical representations (the interested reader could refer to the comprehensive reviews of [START_REF] Oster | Network thermodynamics: dynamic modelling of biophysical systems[END_REF][START_REF] Perelson | Chemical reaction dynamics. part ii: Reaction networks[END_REF]).

In this paper we generalize Mason-Coates type formulas from electrical circuits [START_REF] Mason | Feedback theory -some properties of signal flow graphs[END_REF][START_REF] Coates | Flow-graph solutions of linear algebraic equations[END_REF], providing graphical representations of static susceptivities of gene networks.

We also exploit useful connections between the static response of gene networks and the Dirichlet-to-Neumann map [START_REF] Curtis | Inverse problems for electrical networks[END_REF] which is a relatively new field in graph theory. This field originated from a problem by Calderon [START_REF] Calderon | On an inverse boundary value problem[END_REF]) who asked whether it is possible to reconstruct the conductivity of a conducting plate from measurements of injected currents and voltages on its boundary. This problem gave rise to non-trivial developments in graph theory. Our most important result is related to the Dirichletto-Neumann map and consists in a systematic use of the notion of boundary. Boundaries serve to delimitate subgraphs or to couple subgraphs one with another like in [START_REF] Kholodenko | Untangling the wires: a strategy to trace functional interactions in signaling and gene networks[END_REF]. We define response of subgraphs as solutions of the Dirichlet or the Neumann problems. The novelty of our method is that boundaries delimitating analysed subgraphs are chosen according to the missing piece of information.

This paper is organized as follows. In the first section we present the interaction graph, the subgraphs and their boundaries. In the second section we define linear response and discuss an electrical network analogy. In the third section we construct Mason-Coates solutions to the linear Dirichlet and Neumann problems and we extend the applicability of the results to nonlinear response. Finally, we give two biological applications in the last section.

Steady state shifts, interaction graph

The state of a biological network is given by a vector X ∈ R n , whose components are concentrations of various interacting actors, such as DNA regions coding for genes, RNA transcripts, various produced and regulating proteins, metabolites.

Here we consider a differential dynamics for X, dX dt = F (X, P ), F : R n × R k → R n , where P ∈ R k stems for a set of parameters. This assumption is rather general, because boolean and piecewise-deterministic dynamics can be approximated by differential inclusions [START_REF] Gouzé | A class of piecewise linear differential equations arising in biological models[END_REF]. Although not proven in full generality but used in practice, dynamical systems with delays can be approximated by differential systems either by introducing extra variables [START_REF] Belych | On the problem of replacing an ordinary differential equation with constant time delay by the system of ordinary differential equations[END_REF] or by centre manifold techniques [START_REF] Wang | Dimensional reduction for nonlinear time-delayed systems composed of stiff and soft substructures[END_REF]. The main restriction in our theory is that we consider nondegenerate stable steady states, which are hyperbolic fixed points of the dynamics. Steady states are solutions of the nonlinear system of equations:

F (X, P ) = 0 (1)
Changes in the control parameters P produce steady state shifts. If det[J(X 0 , P 0 )] = 0, where J ij = ∂Fi ∂Xj is the Jacobian calculated at a non-degenerate stable steady state X 0 , then from the implicit function theorem Eq.1 defines X as a function of P locally, in a neighborhood of (X 0 , P 0 ). This local dependence represents the static linear response of the system. In section 3.4 we shall give a global existence theorem for such a function.

Interaction graph

The Jacobian J calculated in a state X introduces in a natural way a signed oriented graph, called the interaction graph which is a triplet (I, E, σ). The set of nodes I = {1, . . . n} consists of the interacting actors. The set of oriented edges A ⊂ I × I is defined by (j, i) ∈ A, iff J ij = 0. The edge sign is a function σ : A → {-1, 1}, σ(j, i) = sign[J ij ]. An edge (j, i) is called positive whenever σ(j, i) > 0 and negative whenever σ(j, i) < 0.

From a biological point of view there is an arc from j to i if the actor j has a direct influence on the dynamics of i. This influence can be positive (activation) or negative (repression). Many examples can be considered. For instance j may be a transcription factor regulating the expression of i, or a protein involved in the phosporylation or methylation of i, or an enzyme involved in the production of i, etc.

When F is non-linear, the signs of its partial derivatives may change, therefore the interaction graph generally depends on the state X where it is calculated. This is a weakness of the interaction graph compared with other more complex graphical representations such as "bond graphs" [START_REF] Perelson | Chemical reaction dynamics. part ii: Reaction networks[END_REF][START_REF] Oster | Network thermodynamics: dynamic modelling of biophysical systems[END_REF]. Nonetheless, this weakness is largely compensated by the fact that many qualitative properties of the dynamics and of steady states depend on topological conditions on the interaction graph [START_REF] Thomas | On the relation between the logical structure of systems and their ability to generate multiple steady states or sustained oscillations[END_REF][START_REF] Soulé | Graphic requirements for multistationarity[END_REF][START_REF] Gouzé | Positive and negative circuits in dynamical systems[END_REF][START_REF] Snoussi | Necessary conditions for multistationarity and stable periodicity[END_REF].

Boundary Boundary is a key concept of any modular approach trying to reduce the overwhelming complexity of large networks to understandable behaviour of simpler subsystems. Furthermore, data in molecular biology can be extensive but it is rarely complete. Our modular approach to data analysis can be used to fill in missing information, to check existing data, or to correct eventual errors.

Because there is no precise criterion to fix the boundary between the actors that we consider and those that we forget, we shall allow a free choice of the boundary. We shall see that instead of becoming a handicap, this freedom of choice becomes a handy tool.

The construction by which we freely choose a boundary in an interaction graph is the following. Let G be the set of nodes that we isolate in the larger set I. The entrance boundary of G, denoted by in G is the set of nodes of G that have incoming arcs from the exterior, i.e. from nodes of I which are not nodes of G. The other nodes of G that are not directly influenced by the exterior are in the interior of G which is denoted G = G \ in G.

Paths, loops, loop partitions

A path c = i ❀ j is a sequence of nodes (i 1 = i, i 2 , . . . , i p = j) such that (i k , i k+1
) is an edge of the graph and such that all the nodes in the sequence are visited just once. Notice that such a path exists in the interaction graph if J i2i = 0, . . . , J jip-1 = 0. A path i ❀ j in the interaction graph is called positive if the path product of the Jacobian elements a ij = p-1 k=1 J i k+1 i k is positive (negative). A loop l = i ❀ i is a sequence of nodes {i 1 = i, i 2 , . . . , i p = i} such that (i k , i k+1 ) is an edge of the graph and such that all the nodes in the sequence are visited just once with the exception of the two terminals that coincide. A loop is called positive (negative) if the product of the Jacobian elements along it is positive (negative). Self-interaction loops (i, i) contain only one node and correspond to diagonal elements of the Jacobian J ii .

A loop partition of a graph is a partition of the nodes into disjoint loops (including self-interaction loops). The set of loop partitions of a subgraph G is denoted L(G). For L ∈ L(G), |L| denotes the number of loops in L. For instance if G = {1, 2, 3} and the oriented edges are

{(1, 1), (2, 2), (3, 3), (1, 2), (2, 3)}, L(G) consists of two loop partitions L 1 = {(1, 1), (2, 2), (3, 3)}, |L 1 | = 3 and L 2 = {(1, 2, 3)}, |L 2 | = 1.
2. Linear response of a network to influences from the exterior

Electrical networks, flowgraphs, Markov chains

Classical examples of linear networks are the linear electrical networks. Let us consider an electrical network with n nodes. Nodes i and j are connected by wires having admittances Y ij = Y ji . Injecting currents I i in some or all of the nodes we produce the steady state potentials V i . Applying Ohm's law and Kirchhoff's first law to the node i it follows that j =i Y ji (V j -V i ) + I i = 0. In matrix form, the linear relation between node voltages and node current sources reads:

Ỹ V = I.
(2)

Ỹ is the node admittance matrix and is obtained from the edge admittances:

Ỹij = -Y ij if i = j j =i Y ij if i = j.
(3)

In order to express the voltages for a given current configuration we need to solve the system of Kirchhoff-Laplace equations 2.

Notice that the matrix Ỹ is singular: (1, . . . , 1) ∈ Ker( Ỹ ). This comes from a special symmetry of electrical networks meaning that voltages are determined up to a constant and that the only measurable quantities are voltage differences.

There are many other examples coming from different fields of science were linear equations are interpreted as graphs.

Markov or semi-Markov processes on multistate stochastic networks occur frequently in biology and physics and can be used for data analysis. For instance the evolution of a disease can be treated as a series of stochastic transitions between various grades of disease [START_REF] Yau | Analysis of censored and incomplete survival data using flowgraph models[END_REF]. The possible transitions and transition probabilities can be gathered on a flowgraph. A flowgraph is a weighted oriented graph. Each edge is labeled by a transition rate (or transmittance, or gain) which is the probability p ji to perform the jump from the state j to the state i divided by the mean waiting time τ ji from a state j to i: t ji = p ji /τ ji . The probabilities π i of being in a state i satisfy the equilibrium equation k π k t ki = ( k t ki )π i which is the same as Eq. 2 with zero currents I = 0.

Analogy between steady state shifts and electrical networks

In order to obtain an analogy between linear response of electrical and biological networks, let us consider that there is a non-degenerated stable state X 0 satisfying Eq.1 for parameters P 0 . Supposing that dynamics is structurally stable the theorem on persistence of hyperbolic sets [START_REF] Smale | Mathematics of Time[END_REF][START_REF] Ruelle | Elements of differentiable dynamics and bifurcation theory[END_REF] implies that small variations of the parameters in Eq.1 produce shifts of this state without loss of stability. Differentiating the equilibrium equations for the nodes in a subgraph G (including the boundary nodes) we obtain:

∀i ∈ G, j∈G ∂F i ∂X j δX j + k∈I\G ∂F i ∂X k δX k + ∂F i ∂P δP = 0.
(4) Let us also consider that all the exterior influences on G are transmitted via its boundary. This means that:

∂F i ∂P = 0, ∀i ∈ G. (5)
The condition in Eq.5 is a modelling assumption, that can be lifted with minor changes. From a statical point of view the parameters P play the same role as the other variables. If the modeler judges that Eq.5 is not reasonable, it is better to consider P as a variable and to connect it to the other variables.

Noticing that for all k ∈ I \ G and i ∈ G, it follows:

∀i ∈ G, j∈G ∂F i ∂X j δX j = -δX f i , (6)
where δX f i denote the "forcing variations" which are non-vanishing only on the boundary, and are defined by:

δX f i = k∈I\G ∂Fi ∂X k δX k if i ∈ in G 0 if i ∈ G. (7)
In matrix form Eq.6 reads:

(-J G )δX = δX f . (8)
Eqs.8 and 2 are analogous. They represent the linear response of a network to influences from the exterior:

• the opposite Jacobian -J G restricted to G and calculated at a non-degenerated stable steady state is analogous to the node admittance matrix Ỹ , • the concentration variations δX are analogous to the voltages V , • the forcing variations δX f are analogous to the injected currents I.

Eq.5 implies that the forcing variations are vanishing in the interior and are non-zero on the boundary of the subgraph G. For electrical networks this means that the influence of the exterior on a subnetwork can be represented by current sources on its boundary.

Electrical networks are different from biological networks because they have much more symmetries, such as:

• The matrix -J G is not necessarily symmetric, while the matrix Ỹ is always symmetric. The interaction graph is intrinsically oriented, while an electrical network is not oriented (a wire conducts in both directions). • The diagonal elements of the matrix Ỹ are obtained from the non-diagonal ones: Ỹii =j =i Ỹij . The diagonal elements of the matrix J G are independent of the non-diagonal ones.

Topological expressions for the linear response of networks

Throughout this section we consider that the Jacobian J G is calculated at a non-degenerate steady state (stability can be imposed but it is not compulsory). For subgraphs, steady states are constrained by imposed boundary variations or forcings. To simplify notation we drop the index G in the Jacobian; it will be supposed that the J is restricted to a subset G, unless otherwise stated. Supposing that J is invertible, solutions of Eqs.8 can be obtained from the inverse matrix J -1 . In this section, the elements of J -1 are connected to the topology of the interaction graph.

Moduli and loops

We start with the following well known relation [START_REF] Bloom | Linear algebra and geometry[END_REF]:

J -1 ij = (-1) i+j ∆ ji ∆ (9)
where ∆ ji is the minor obtained by deleting row j and column i in J, and ∆ = det(J).

Next, we develop the minor ∆ ji into a sum of principal minors [START_REF] Bloom | Linear algebra and geometry[END_REF]:

∆ ji = (-1) i+j j❀i (-1) lj❀i a j❀i ∆ j❀i if i = j ∆ j if i = j.
(10) where j ❀ i is any path in the subgraph G leading from j to i, l j❀i is the number of edges in the path, a j❀i is the product of elements of J along this path, ∆ j❀i is the principal minor obtained by deleting all the rows and columns whose indices are included in the path, ∆ j is the principal minor obtained by deleting the line and the column j. Finally principal minors of the Jacobian can be related to loops in the subgraph G:

∆ j❀i = L∈L(Gj❀i) (-1) dim(∆j❀i)-|L| lp(L) (11) ∆ i = L∈L(Gi) (-1) dim(∆i)-|L| lp(L) (12) ∆ = L∈L(G) (-1) dim(∆)-|L| lp(L). ( 13
)
where lp(L) is the loop product defined as the product of elements of J along the loops of L. G j❀i is the subgraph obtained by deleting from G all nodes in the path j ❀ i, while G i is obtained by deleting node i. dim denotes the dimension of the minor. We now introduce a quantity later referred to as modulus which measures rigidity of the network and can be seen as the inverse of sensibility, as shown later:

Definition 1. Let G be a subgraph of I, let ∆, ∆ i , ∆ j❀i be minors of the Jacobian restricted to G and defined as above. The modulus of a path j ❀ i is:

C j❀i = (-1) lj❀i+1 ∆ ∆ j❀i = = L∈L(G) (-1) |L| lp(L) L∈L(Gj❀i) (-1) |L| lp(L) . ( 14
)
The modulus of a node i is:

C i = - ∆ ∆ i = L∈L(G) (-1) |L| lp(L) L∈L(Gi) (-1) |L| lp(L) . ( 15
)
We shall usually consider non-degenerate steady states and ∆ = 0, therefore moduli never vanish; they can instead diverge when minors of J vanish.

Neumann and Dirichlet problems

In electrical networks, fulfillment of Kirchhoff-Laplace Eqs. 2 implies a unique set of interior voltages (up to a constant) due to imposed boundary voltages (Dirichlet problem) or a unique set of voltages (up to a constant) everywhere due to imposed boundary currents (Neumann problem). For biological networks, the Neumann and the Dirichlet problems have the following equivalents:

• Linear Dirichlet problem: determine the variations δX i , i ∈ G of the interior nodes, when the values of the variations are imposed (or known) on the boundary δX i , i ∈ in G. For biological networks the Dirichlet problem shows its utility when one possesses only a partial knowledge of the system [START_REF] Kuipers | Qualitative reasoning. Mdeling and simulation with incomplete knowledge[END_REF][START_REF] Siegel | Qualitative analysis of the relation between dna microarray data and behavioral models of regulation networks[END_REF]. We can define the corresponding Dirichlet static susceptivitities:

χ D ij = ∂X i ∂X j , i ∈ G, j ∈ in G. (16)
• Linear Neumann problem: determine the variations δX i , i ∈ G everywhere, when the forced variations are imposed on the boundary δX f i , i ∈ in G. Unlike injected currents forced variations (defined by Eqs.7) can not be measured directly. In electrical networks currents can be measured as voltages across calibrated resistors of apparatuses coupled to the network. We can not see what can be a calibrated molecular apparatus for biological networks. In spite of this drawback, the Neumann problem is important in theory of control of biological networks: it describes shifts of equilibria under constraints imposed by the change of external conditions. The corresponding static susceptivities are defined as:

χ N ij = ∂X i ∂X f j , i ∈ G, j ∈ in G. (17)
The following theorem states the existence and gives a solution to the Neumann problem.

Theorem 1. Let G be a subgraph of I. Let J be the restriction of the Jacobian of F to G, calculated at a non-degenerated, Neumann constrained (given forcings) steady state: J ij = ∂Fi ∂Xj , i, j ∈ G. Let P G denote the set of paths included in G. If det(J) = 0 and if there is no direct influence of the parameters on the nodes of G, i.e. ∂Xi ∂P k = 0, ∀i ∈ G, then the response of i ∈ G to small changes of the exterior of G satisfies:

if i ∈ G, δX i = j∈ in G j❀i∈PG a j❀i C j❀i δX f j , ( 18
) if i ∈ in G, δX i = δX f i C i + j ∈ in G, j = i j❀i∈P G a j❀i C j❀i δX f j . ( 19
)
Remark 1. We have supposed that ∆ = det(J) = 0. According to Eqs.14,15, we can supplement Eqs.18,19 with the following formal rule coping with diverging moduli: 1/∞ = 0.

Remark 2. Using a language coming half from mechanics, half from flowgraphs we can say that a "force" a j❀i δX f j propagates from the boundary node j along the path j ❀ i. This force is bigger when the product of interaction coefficients ∂Fi k+1 ∂Xi k along the path is bigger. C j❀i is the ratio force/response and therefore can be called "path modulus ". A large path modulus implies a small response at the end of the path, even if the force is big. Therefore, the modulus is the inverse of sensitivity.

Eq.18 expresses interior δX i as functions of the boundary forcings δX f j . Eq.19 gives the relation between forcings and variations on the boundary. The corresponding linear mapping Λ

N D : R m → R m (m = # in G, Λ N D ij = j❀i∈PG aj❀i Cj❀i + 1 Ci δ ij )
is the equivalent of the Neumann-to-Dirichlet map for electrical networks [START_REF] Curtis | The dirichlet to neumann map for a resistor network[END_REF]. The inverse (if it exists) of this mapping is the Dirichlet-to-Neumann map [START_REF] Curtis | The dirichlet to neumann map for a resistor network[END_REF].

Theorem 1 implies the following:

Corollary 1. If there are no paths in P G connecting two different nodes on the boundary (boundary is not necessarily disconnected because connections passing through the exterior are allowed), and if for all nodes i ∈ in G, moduli C i are finite and non-vanishing, then the Neumann-to-Dirichlet mapping is diagonal and its inverse exists. In this case, the forcings are proportional to the variations on the boundary:

δX f i = C i δX i , i ∈ in G (20)
By using the Dirichlet-to-Neumann map we can associate to any solution of the Neumann problem, a solution of the Dirichlet problem. We can thus obtain formulae for the solutions of the Dirichlet problem. In [START_REF] Siegel | Qualitative analysis of the relation between dna microarray data and behavioral models of regulation networks[END_REF] we have used a direct method to obtain solutions to the Dirichlet problem. Briefly, we have used equilibrium equations for the interior nodes in the same way as in Theorem 1 (which uses equilibrium equations for all the nodes including the boundary). Let us state the result:

Theorem 2. Let G be a subgraph of I. Let G = G \ in G be the interior of G. Let J be the restriction of the Jacobian of F to G calculated in a non-degenerated, Dirichlet constrained (given boundary values, free interior) steady state: Ji,j = ∂Fi ∂Xj , i, j ∈ G. Let us suppose that ∆ = det( J) = 0 and that there is no direct influence of the parameters on the interior nodes of G, i.e. ∂Fi ∂P k = 0, ∀i ∈ G.

Then, the response of i ∈ G to small changes on the boundary of G is given by:

δX i = j∈ in G j❀i∈P G a j❀i Cj❀i δX j . (21)
• P G denotes the set of paths included in G, starting on the boundary and that do not return to the boundary.

• Cj❀i = (-1) l k(j)❀i +1 ∆ ∆k(j)❀i
denotes the path modulus of the internal path k(j) ❀ i where k(j) ∈ G is the second node after j of the path j ❀ i. If k(j) = i, then Cj❀i = Ci = -∆ ∆i Remark 3. Theorems 1,2 make no assumption on the connectivity of the subgraph G. For disconnected internal nodes we should have δX i = 0, meaning that at steady states X i does not depend on the parameters. For such a node the stationarity equations read F i (X i ) = 0.

Eqs.14,15 give the loop products expression of moduli. They generalize Mason and Coates gain formulae [START_REF] Mason | Feedback theory -some properties of signal flow graphs[END_REF][START_REF] Coates | Flow-graph solutions of linear algebraic equations[END_REF] from electrical networks. Susceptivities are related to moduli according to the following equations:

χ N ij = j❀i∈PG a j❀i C j❀i + 1 C i δ ij (22) χ D ij = j❀i∈P G a j❀i Cj❀i ( 23
)
Let us notice that the restriction of the Neumann susceptibility matrix to the boundary ((χ

N | in G ) ij = χ N ij , i, j ∈ in G)
is the Neumann-to-Dirichlet map.

Signs of moduli

The signs of moduli are important for qualitative discussions of the transmission of influences [START_REF] Siegel | Qualitative analysis of the relation between dna microarray data and behavioral models of regulation networks[END_REF]. Eqs.22,23 lead to an interesting possibility. Even if all the paths from j to i correspond to globally positive regulation (path products a j❀i are all positive), it is possible to have negative susceptivity χ ij , which means negative correlation between X i and X j . This possibility occurs when the moduli C j❀i are negative (see Section 4.2) and is a feed-back effect. Indeed, from Eqs.14,15 it follows:

Property 1. If the subgraph G contains no positive loops for a given non-degenerate constrained steady state, then all its path moduli C j❀i and node moduli C i are positive for that state.

Significance and limits of the approach

Eqs. 18,19,21 should be interpreted as differential constraints connecting small variations of the concentrations at steady state, when this steady state is shifted as a result of the change of parameters. The significance of these constraints is the propagation of influence along the interaction graph: for any analyzed subset of the interacting species concentration variations of internal species are obtained from concentration variations of boundary species. Mathematically, these constraints express a local functional relation between the internal variables and the boundary variables, which follows from the implicit function theorem. Let us notice that we have not addressed yet the question of the global existence of an unique implicit function expressing the internal variables as functions of the boundary variables (the non-linear Dirichlet problem). The existence of such a function depends on the non-linearities of the system and can not be proven without further additional hypothesis. Without being exhaustive, let us give here some results applying to particular cases.

Let us consider that for certain subsets G ⊂ I we are able to specify global interiors, i.e. sets of species that are never directly influenced by external species, whatever the concentrations. The steady state equations for internal species constrained by boundary species read:

F i (X, X ′ ) = 0, i ∈ G (24)
where X,X ′ have coordinates in G, in G, respectively. The implicit function theorem implies that X is a function of X ′ in a neighbourhood of any nondegenerate steady state satisfying det( J) = 0. The global existence of the implicit function follows from an existence result for the solution of the system 24 and from the following global univalence theorem ensuring the uniqueness of this solution:

Theorem 3. Let us consider that F i (X, X ′ ) = Φ i (X, X ′ ) -λ i X i , i ∈ G where λ i > 0 and Φ i (X, X ′ )
are differentiable, bounded and satisfy Φ i (. . . , X i = 0, . . . , X ′ ) > 0 (25)

Then for any X ′ the system 24 has at least a solution X such that all concentrations X i , i ∈ G are positive.

Let us consider that the interaction subgraph corresponding to the interior of G has no positive loop, whatever the concentrations may be. Then the solution of the system 24 is unique for any X ′ .

We shall only give a brief outline of the proof. The first part of the Theorem 3 is based on the following standard mathematical lemma which is a consequence of Poincaré-Hopf formula:

Lemma 1. Let D be a smooth ball in a metric space, let S be the boundary of D. Let F be a vector field defined on a neighbourhood of D. If F points inward D at any point of S, then F admits a zero in the interior of D.

The second part of the theorem is known as the Thomas's conjecture [START_REF] Thomas | On the relation between the logical structure of systems and their ability to generate multiple steady states or sustained oscillations[END_REF][START_REF] Gouzé | Positive and negative circuits in dynamical systems[END_REF][START_REF] Snoussi | Necessary conditions for multistationarity and stable periodicity[END_REF][START_REF] Soulé | Graphic requirements for multistationarity[END_REF] and can be proven as follows. When the subgraph G has no positive loops whatever the concentrations, Eqs.11,12,13 imply that all the minors of the opposite Jacobian -J are positive. From this and from the Gale-Nikaido-Inada theorem [START_REF] Parthasarathy | On Global Univalence Theorems[END_REF]it follows the uniqueness of the equilibrium. The complete proof of Thomas's conjecture can be found in [START_REF] Soulé | Graphic requirements for multistationarity[END_REF].

The conditions of the first part of the Theorem 3 are fulfilled by gene networks. Indeed, the terms λ i X i correspond to degradation or to dilution produced by cell growth. The production terms Φ i are bounded because they saturate at large concentrations. Eq.25 is natural, meaning that concentrations can never become negative.

We come now to the experimental significance of our approach. The type of experiment that we consider is the following: the set of parameters is changed such that we start and we end in a non-degenerate stable steady state. If the conditions of Theorem 3 are fulfilled then the nonlinear Dirichlet problem has a solution and the internal variables are constrained to be functions of the boundary variables. The finite variations of these functions between the start and the end of the experiment can be obtained by integrating the differential form in Eq.21 along any differentiable path D connecting the two steady states and made of nondegenerate constrained steady states (non-degenerate Figure 1. Cusp catastrophe. Two trajectories T1 and T2 start in the steady state P and end in the steady states P1 and P2 respectively. In order to find the variation of the internal variable X 3 one could integrate the differential response along the smooth branch P P 1 (for T 1 ) and along the branches P P 1 P ′ 1 and P ′ 2 P 2 (for T 2 ). For T 2 the finite jump P ′ 1 P ′ 2 should be added.

solutions of Eq.24)

.

∆X i = D j∈ in G j❀i∈P G a j❀i Cj❀i dX j . ( 26
)
The uniqueness of the implicit function ensures the path independence of the result. We must emphasize that Eq.26 is not destined to provide numerical values. It serves to identify the signs of different contributions to the concentration variations and leads to qualitative equations in sign algebra [START_REF] Kuipers | Qualitative reasoning. Mdeling and simulation with incomplete knowledge[END_REF] (details can be found in [START_REF] Siegel | Qualitative analysis of the relation between dna microarray data and behavioral models of regulation networks[END_REF]). Similar results hold for shifts of Neumann constrained steady states (nonlinear Neumann problem).

This result fails when the steady state shift meets a singular point, defined by vanishing moduli (infinite sensitivity). Then the steady state loses stability. The new attractor is either a distance apart (for instance in the cusp catastrophe), or it has a different type (for instance in the Hopf bifurcation a point attractor becomes a stable limit cycle) from the old one.

The cusp catastrophe occurs rather often in molecular biology. It occurs in the well studied lactose operon, that we shall present in the next section. This situation is represented in Fig. 3.4. The boundary variables X ′ 1 , X ′ 2 influence the internal variable X 3 . Suppose that X ′ 1 changes. Then the state of the system moves along the trajectories starting from P . For small changes ending in P 1 there is a differentiable path of stable steady states connecting P and P 1 and we can apply Eq.26. For large changes ending in P 2 there is no such path. The steady state shift is discontinuous: the state jumps from the attractor branch P P 1 P ′ 1 to the branch P 2 "P ′ 2 P 2 . Eq.26 has to be changed by adding a finite jump X 3 (P ′ 2 ) -X 3 (P ′ 1 ).

Examples

Lactose operon

The main enzymes for the lactose (L) metabolism in E.coli are LacY (lactose permease) allowing the uptake of lactose, LacZ (β-galactosidase) catalyzing the degradation of lactose to glucose [START_REF] Mackey | Modeling operon dynamics : the tryptophan and lactose operons as paradigms[END_REF][START_REF] Yildirim | Feedback regulation in the lactose operon: a mathematical modeling study and comparison with experimental data[END_REF]. The transcriptional regulators for the genes lacY and lacZ are an activator (CRP) and a repressor (LacI). An inducer (cAMP) binds to the activator and triggers it. Lactose, under an isomeric form named allolactose, binds to the repressor and inhibits it. The glucose inhibits the activator. The interaction graph for this system is represented in Fig. 2. The exterior and interior lactose are denoted L e , L i , respectively. Negative self-interaction loops should be added to each node, in order to take into account degradation processes or dilution produced by cell growth. These correspond to negative diagonal elements for the Jacobian (self-interactions).

There are four loops in the interaction graph: two positive (even number of negative regulations) ones l + 1 = {L i , LacI, LacY, L i }, l + 2 = {L i ,G,cAM P , CRP, LacZ, L i }, and three negative (odd number of negative regulations) ones l - 1 = {L i , LacI, LacZ, L i }, l - 2 = {G, cAM P, CRP, LacZ, G}, l - 3 = {L i , G, cAM P, CRP, LacY, L i }. The existence of positive loops is a necessary condition [START_REF] Soulé | Graphic requirements for multistationarity[END_REF] for the observed bistability of the operon lactose.

Experiments show that when L e is increased, the operon switches from a transcriptionally blocked, lactose poor state to a transcriptionally active, lactose rich state. We intend to check whether the linear response of the model given in Fig. 2 is coherent to observed variations, namely whether the Dirichlet susceptivity

χ D L i ,L e is positive. From Theorem 2, χ D L i ,L e = a L e ❀L i CL e ❀L i
, where a L e ❀L i > 0 from the construction of the model (see Fig. 2). This means that the influence of L e on L i is positive provided that the modulus of the path L e ❀ L i is positive. The modulus is given by:

CL e ❀L i = = χ L i 1 - 2 i=1 lp(l + i ) - 3 i=1 lp(l - i ) + lp(l + 1 ) lp(l - 2 ) 1 -lp(l - 2 ) (27)
where χ L i is the absolute value of the self-interaction on L i , lp is the loop product divided by the product of absolute values of self-interactions for nodes in the Static response of interaction networks O.Radulescu et al loop. lp is positive for positive loops and negative for negative loops.

The condition for a monotonic response is given by the following: Property 2. The dependence of L i on L e is monotonically increasing provided that the modulus CL e ❀L i is positive, i.e.:

lp(l + 1 )[1 -lp(l - 2 )] + lp(l + 2 ) < 1 - 3 i=1 lp(l - i ). (28)
The monotonicity condition 28 is valid for small positive loops products. The typical scenario for strong positive loops is the cusp catastrophe: the operon becomes bistable and its response is discontinuous and hysteretic [START_REF] Mackey | Modeling operon dynamics : the tryptophan and lactose operons as paradigms[END_REF]. Bistability requires the existence of a branch of states with negative moduli on which the condition 28 is broken, but this branch is usually unstable (see also next subsection).

The susceptivity approach can predict the effect of gene knock-outs. For instance, a knock-out of the inducer or of the activator cancels the loops l + 2 , l - 2 , l - 3 . According to Eq.27, the modulus increases, therefore the uptake of lactose decreases, provided that lp(l + 2 ) > -lp(l - 2 ) lp(l - 1 ) -lp(l - 3 ). Similarly, a knock-out of the repressor cancels the loops l + 1 , l - 1 . The uptake of lactose decreases provided that lp(l + 1 ) > lp(l + 1 ) lp(l - 2 ) -lp(l - 1 ).

Negative moduli and non-monotonic response

For the lactose operon the input/output relation relating outside and inside lactose is monotonic and the moduli for stable branches are positive. We build here an artificial example of switch with non-monotonic response.

Let us suppose that X 1 regulates positively X 2 and that the modulus C12 < 0. Then increasing X 1 produces a decrease in X 2 . In order to create an example of this kind one should consider at least another component X 3 of the system. Indeed, stability asks for eigenvalues of the Jacobian to be in the left half of the complex plane; in dimension 2 this is equivalent to the absence of positive loops, hence moduli are positive at steady states (Property 1). A simple example illustrating the possibility of having negative moduli is:

   dx1 dt = -x 1 + 1 + λ, dx2 dt = -2x 2 -x 3 + x 1 , dx3 dt = f (x 3 ) + 3x 2 (29) with f (x) = -ax(x -1)(x -2).
The steady state is shifted by changing the parameter λ, which at equilibrium is equivalent to changing the boundary value

x 1 = 1 + λ.
The solution of the Dirichlet problem is δx 2 = δx1

C12

,

δx 3 = δx1

C13

, where

C12 = 3 -2 df dx 3 / - df dx 3 (30) C13 = 3 -2 df dx 3 . ( 31
) 0 1 2 -0.5 0 0.5 1 1.5 2 C 12 < 0 x 2 x 3 x 1 x 2 ,x 3 C 12 < 0 0 1 2 -0.5 0 0.5 1 1.5 2 C 13 < 0 x 2 x 3 x 1 x 2 ,x 3 C 13 < 0 a) b)
Figure 3. Solution of the Dirichlet problem for the model in Eq.4.2: stable equilibrium values of x 2 and x 3 as a function of the boundary value x 1 . a) Monostable (a = 0.5) with a region (between dotted vertical lines) of negative modulus C12 < 0. b) Bistable (a = 2) with discontinuous, hysteretic response; the response of x2 is non-monotonic, the response of x3 is monotonic. C13 can be negative when x3 is between the two dotted horizontal line, but the corresponding states are unstable.

The model has unique equilibrium for a < 3/2 and is bistable for a > 3/2. In the monostable regime (a < 3/2) C13 > 0 everywhere; C12 ≤ 0 when x 3 is inside the interval

I 1 = [1 -1/ √ 3, 1 + 1/ √ 3] and is positive outside the interval I 1 . For x 3 = 1 -1/ √ 3, 1 + 1/ √ 3, 1 C12
and therefore the susceptivity χ D 21 vanishes: the response curve for x 2 represented in Fig. 4.2 have a maximum and a minimum at these points.

For a > 3/2, the system is bistable. We have C12 > 0 unless when

x 3 ∈ I 2 = [1 -1/ √ 3, 1 - 1/3 -1/(2a)] ∪ [1 + 1/3 -1/(2a), 1 + 1/ √ 3] and C13 > 0 unless when x 3 ∈ I 3 = [1 -1/3 -1/(2a), 1 + 1/3 -1/(2a)].
The region I 3 where C13 < 0 is in fact unstable and is never reached. When x 1 increases x 3 jumps discontinuously from one attractor to another without taking values inside the interval I 3 . The jump occurs at different values of x 1 on increase and decrease of the control parameter. Thus the response of x 3 is discontinuous and hysteretic, but monotonic. The response of x 2 is also discontinuous and hysteretic, but non-monotonic. There is a small region on the response curve of x 2 where C12 < 0; this region corresponds to the intersection between the stability domain for x 3 and I 2 .

Regulation of lipogenesis in hepatocytes

Interaction model Two ways of production of fatty acids coexist in liver. Saturated and mono-unsaturated fatty acids (MUFA) are produced from citrates thanks to a metabolic pathway composed of four enzymes, namely ACL (ATP citrate lyase), ACC (acetyl-Coenzyme A carboxylase), FAS (fatty acid synthase) and SCD1 (Stearoyl-CoA desaturase 1). Polyunsaturated fatty acids (PUFA) are synthesized from essential fatty acids provided by nutrition; D5D (Delta-5 Desaturase) and D6D (Delta-6 Desaturase) catalyze the key steps of the synthesis of PUFA.

PUFA plays pivotal roles in many biological functions; among them, they regulate the expression of genes that impact on lipid, carbohydrate, and protein metabolism. The effects of PUFA are mediated either directly through their specific binding to various nuclear receptors (PPARα -peroxisome proliferator activated receptors, LXRα -Liver-X-Receptor α, HNF-4α) leading to changes in the trans-activating activity of these transcription factors; or indirectly as the result of changes in the abundance of regulatory transcription factors (SREBP-1c -sterol regulatory element bindingprotein-, ChREBP, etc.) [START_REF] Jump | Fatty acid regulation of gene transcription[END_REF].

We consider in our model nuclear receptors PPARα, LXRα, SREBP-1c (denoted by PPAR, LXR, SREBP respectively in the model), as they are synthesized from the corresponding genes and the trans-activating active forms of these transcription factors, that is, LXR-a (denoting a complex LXRα:RXRα), PPAR-a (denoting a complex PPARα:RXRα) and SREBP-a (denoting the cleaved form of SREBP-1c). SCAP -(SREBP cleavage activating protein) is a key enzyme involved in the cleavage of SREBP-1c. We also include in the model "final" products, that is, enzymes ACL, ACC, FAS, SCD1 (implied in the fatty acid synthesis from citrate), D5D, D6D (implied in PUFA synthesis) as well as PUFA themselves.

Relations between the variables are the following. SREBP-a is an activator of the transcription of ACL, ACC, FAS, SCD1, D5D and D6D [START_REF] Jump | Fatty acid regulation of gene transcription[END_REF]. LXR-a is a direct activator of the transcription of SREBP and FAS, it also indirectly activates ACL, ACC and SCD1 [START_REF] Steffensen | Putative metabolic effects of the liver x receptor (lxr)[END_REF]. Notice that these indirect actions are kept in the model because we don't know whether they are only SREBP-mediated. PUFA activates the formation of PPAR-a from PPAR, and inhibits the formation of LXR-a from LXR as well as the formation of SREBP-a (by inducing the degradation of mRNA and inhibiting the cleavage) [START_REF] Jump | Fatty acid regulation of gene transcription[END_REF]. SCAP represents the activators of the formation of SREBP-a from SREBP, and is inhibited by PUFA. PPAR directly activates the production of SCD1, D5D, D6D [START_REF] Miller | Peroxisome proliferators induce mouse liver stearoyl-coa desaturase 1 gene expression[END_REF][START_REF] Tang | Regulation of human delta-6 desaturase gene transcription: identification of a functional direct repeat-1 element[END_REF][START_REF] Matsuzaka | Dual regulation of mouse delta(5)-and delta(6)-desaturase gene expression by srebp-1 and pparalpha[END_REF].

The interaction graph for this model is shown in Fig. 4. Like in the lactose operon example, for each node we have supposed the existence of negative selfinteraction loops.

Fasting-refeeding protocols The fasting-refeeding protocols are suited for studying lipogenesis regulation; during an experimentation, animals (rodents or chicken) were kept fasting during several hours and then refed. Hepatic mRNA of LXR, SREBP, PPAR, ACL, ACC and SCD1 were quantified by DNA microarray analysis. PUFA variations were determined by biochemical measurements.

A compilation of recent literature on lipogenesis regulation provides results of such protocols: SREBP, ACL, ACC, FAS and SCD1 decline in liver during fasting [START_REF] Liang | Diminished hepatic response to fasting/refeeding and liver x receptor agonists in mice with selective deficiency of sterol regulatory elementbinding protein-1c[END_REF]; this state is characterized by an inhibition of fatty acid synthesis and an activation of the fatty acid oxydation. However, [START_REF] Tobin | Cross-talk between fatty acid and cholesterol metabolism mediated by liver x receptor-alpha[END_REF]) showed that fasting rats for 24h increased the hepatic LXR mRNA and Matsuzaka et al [START_REF] Matsuzaka | Dual regulation of mouse delta(5)-and delta(6)-desaturase gene expression by srebp-1 and pparalpha[END_REF] observe no difference in either the hepatic D5D or D6D mRNA level between fasting and refeeding in normal mouse livers.

PPAR [+] PPAR-a PUFA [+] ACL [-] LXR-a SREBP-a SCAP D5D [0] D6D [0] ACC [-] FAS [-] SCD1 [-] LXR [+] EXT SREBP [-]

Static response and topology

One of the advantages of the approach that we present here is that we have no difficulty in focusing on subgraphs of a large network. In order to illustrate this possibility we considered the following biological questions:

Question 1 Lee et al. studied the fatty acids profiles in triglycerides TG, which are the predominant (> 50%) hepatic lipids and also in phospholipids PL which go into cellular membranes. Moreover PL contribute much less than TG to the total lipid mass. These authors [START_REF] Lee | Requirement of pparalpha in maintaining phospholipid and triacylglycerol homeostasis during energy deprivation[END_REF] show that after 72h of fasting [START_REF] Lee | Requirement of pparalpha in maintaining phospholipid and triacylglycerol homeostasis during energy deprivation[END_REF]) fatty acids profiles do not change significantly in PL, but there is a strong increase of TG and of their fatty acids constituents, in particular PUFA. Let us recall that TG are transient storage forms of fatty acids. Based on these experimental findings, we make the hypothesis of a mass increase of regulating PUFA in the hepatic cell. Can we prove this hypothesis indirectly by using the incomplete model, that is, considering that all variations excluding PUFA are known?

Question 2 The dual regulation of desaturases SCD1, D5D and D6D by SREBP and PPAR is paradoxical because SREBP transactivates genes for fatty acid synthesis in liver, while PPAR induces enzymes for fatty acid oxidation. Furthermore, all three desaturases have similar regulation. Nonetheless, on fasting SCD1 decreases, while D5D, D6D have non significative variations. How does the model cope with this?

In order to answer to the Question 1, let us consider the Dirichlet problem for the subgraph G 1 = { PUFA,LXR, LXR-a, SREBP }. From Theorem 2 and the paths labels from Fig. 4 it follows:

δSREBP = ãp δPUFA + ãp ′ δLXR. (32)
Because ãp < 0, ãp ′ > 0, δSREBP < 0, δLXR > 0, the only possible variation compatible with Eq.32 is δPUFA > 0. Applying a similar reasoning to the subgraph G 2 = { PUFA, LXR, LXR-a, SREBP, SREBP-a, SCAP, FAS } we have: Property 3. The only possible variation compatible with δSREBP < 0 and δLXR < 0 is δPUFA > 0. The only possible variation compatible with δFAS < 0 and δLXR < 0, is δPUFA > 0. This proves our hypothesis that regulating PUFA increase during fasting.

In order to answer to the Question 2, let us see how the model connects the variations of PUFA and D6D. PUFA being on the boundary we need to consider the Neumann problem. Let us consider the subgraph G 3 = {X=PUFA, W =PPAR, PPAR-a, Z=D6D, SREBP-a, SREBP, LXR-a, U =LXR, SCAP}. Information from litterature suggests that G 3 has the boundary in G 3 = {U, W, X}, meaning that all external interaction on G 3 acts via LXR, PPAR, PUFA. Nevertheless, in order to understand the differential behaviour of D6D and SCD1, let us suppose that D6D has unknown extra regulation that is not represented. This can be taken into account easily by including D6D in the boundary:

in G 3 = {U, W, X, Z}.

From Theorem 1 and the paths labels from Fig. 4 it follows:

δX = 1 C [ δX f χ X + ãp3 δZ f + ãp9 δW f + ãp10 δU f ] δZ = 1 C [(ã p1 + ãp2 + ãp4 + ãp7 + ãp8 )δX f + + ãp5 δU f + ãp6 δW f + δZ f χ Z ]. ( 33
)
where ã are path products divided by the products of the absolute values of the self-interactions χ of nodes in the paths. The modulus is: C = 1 -lp(l + ) -lp(l - 1 ) -lp(l - 2 ) -lp(l - 3 ) -lp(l - 4 ).

(34) Eqs.33, 34 reveal the importance of the dual regulation of D6D by PPAR and SREBP. The existence of the positive loop l + involving PUFA, D6D and PPAR-a decreases the modulus C in Eq.34 therefore, according to Eq. 33, it boosts PUFA variations. The biological significance of this has been pointed out by Nakamura and Nara [START_REF] Nakamura | Essential fatty acid synthesis and its regulation[END_REF] who refer to this phenomenon as a compensatory reaction to the increased demand of PUFA caused by oxidation at fasting.

In Eq.33 the D6D variation δZ is the sum of different contributions of different signs: the influences transmitted along the paths p 1 , p 2 , p 4 , p 8 are negative and those transmitted along the paths p 7 , p 5 , p 6 are positive. The model is therefore compatible with compensated variations δZ = 0, even in absence of unknown interaction whose effect is the forcing δZ f . Nevertheless, it is difficult to explain the difference between D6D and SCD1 without the extra forcing, because this forcing is the only difference between the two products in the model. [START_REF] Nakamura | Essential fatty acid synthesis and its regulation[END_REF] arrive at a similar conclusion in their analysis of the regulation of D6D. To our knowledge, the nature of this extra regulation is still unknown.

Discussion

Static response provides quantitative and qualitative information on how the steady state is changed by forcings.

We have shown here how to relate static response of a biological network to the topology of its interaction graph. The key ingredient is the use of moduli that express the rigidity (opposed to sensitivity) of the network. Mason-Coates formulae from electrical networks computes moduli from loop products and were applied to biological networks. A general consequence of these formulae is that positive (negative) loops decrease (increase) moduli and rigidity. Static susceptivity was related to moduli by using path expansions. An interesting possibility arises when moduli and susceptivities have opposite signs: in this case the sign of the correlation between two nodes can be opposite to the sign of the regulation paths connecting them.

Susceptivities characterize linear response which generally applies to small variations of the concentrations. We showed that in the absence of bifurcations, local linear response can be integrated in order to obtain nonlinear response. This argument justifies qualitative analysis techniques applying to situations when the interactions have constant signs and when the constrained steady states of analysed subgraphs do not bifurcate under forcings. Our formalism is well adapted for qualitative analysis of differential transcriptional, metabolic or proteomic data. It allows several types of analysis: checking the consistency between model and data, or between different types of data, filling in missing information, and predicting the effects of gene knock-outs or of genetic deficiencies. Qualitative analysis of biological networks has been described in detail in [START_REF] Siegel | Qualitative analysis of the relation between dna microarray data and behavioral models of regulation networks[END_REF].

An important ingredient in our approach is the systematic use of boundaries. This allows to separate subgraphs for independent study. The type of biological problem dictates the choice of the subgraph and the type of boundary conditions (Neumann or Dirichlet).

Boundaries could be useful in further developments such as hierarchical connection of subgraphs. Subgraphs can be connected in series (the entrance boundary of one is included in the exit boundary of another) or in parallel (boundary subsets are common) and their equivalent response can be obtained from Eqs. 18,19. Another promising direction is the reduction of complexity of biological networks. One would like to classify all the graph contractions that preserve the response of the system. This has been done for planar electrical networks [START_REF] Curtis | Circular planar graphs and resistor networks[END_REF] in relation to the Dirichlet-to-Neumann map. There is no similar work for biological networks.

Figure 2 .

 2 Figure 2. Operon lactose interaction network. L int and L ext stand for L i and L e respectively.

Figure 4 .

 4 Figure 4. Interaction graph for a model of regulation of the synthesis of fatty acids. The node EXT represents the exterior world. Negative self-interaction loops on nodes are omitted for sake of clarity. Labelled loops: l - 1 = { PUFA, SCAP, SREBPa, D6D, PUFA}, l - 2 = {PUFA, SREBPa, D6D,PUFA}, l - 3 = {PUFA, LXRa, SREBP, SREBPa, D6D, PUFA}, l - 4 = {PUFA, D6D, PUFA}, l + = { PUFA, PPARa, D6D, PUFA}. Labelled paths p = {PUFA, LXR-a, SREBP}, p ′ = {LXR, LXR-a, SREBP}, p1 = {PUFA, D6D}, p2 = {PUFA, SCAP, SREBPa, D6D}, p3 = {D6D, PUFA}, p4 = {PUFA, LXRa, SREBP, SREBPa, D6D}, p5 = {LXR, LXRa, SREBP, SREBPa, D6D}, p 6 = {PPAR, PPARa, D6D}, p7 = {PUFA, PPARa, D6D}, p8 = {PUFA, SREBPa, D6D}, p9 = {PPAR, PPARa, D6D, PUFA}, p10 = {LXR, LXRa, SREBP, SREBPa, D6D, PUFA}.

Prepared using rsifpublic.cls

Acknowledgements

This research was supported by ACI IMPBio, a French Ministry for Research program on interdisciplinarity. We thank M. Crouzeix, A. Gorban, C. Soulé and J. Scott for inspiring discussions. We thank E. Pécou for suggesting us to use Lemma 1, and the reviewers for constructive comments.