J.-P Morillon

S Doubabi

R Costalat

Nonlinear oscillations of the Goodwin metabolic pathway with coupling

Keywords:

An analytic method has been applied to a biochemical pathway with coupling. The harmonic balancing technique can be used to determine the effects of coupling on the Goodwin metabolic pathway. Its behaviour is modelled by a set of coupled ordinary nonlinear differential equations. This one is transformed to a single feedback loop control with a nonlinear function and a linear filter. It is shown that coupling modifies the linear filter such that the frequency and the amplitude of the oscillations decrease. The results obtained here are derived by both analytical and numerical studies, and compared with previous stability analysis of associated pathways.

Introduction

The most commonly encountered form of regulated biochemical pathway, generally referred to as the Yates-Pardee or Goodwin metabolic pathway, consists of a single pathway of enzymatic reactions, where the last product inhibits the first enzyme (single loop negative feedback) [START_REF] Walter | Kinetic and thermodynamics aspects of biological and biochemical control mechanisms[END_REF][START_REF] Goodwin | Analytical Physiology of Cells and Developing Organisms[END_REF]. When the length of a Goodwin metabolic pathway is increased, the stability domain of its unique steady state is decreased [START_REF] Rapp | Analysis of biochemical phase shift oscillators by a harmonic balancing technique[END_REF][START_REF] Walter | Kinetic and thermodynamics aspects of biological and biochemical control mechanisms[END_REF].

From a different point of view, G. Chauvet has suggested that the association of metabolic pathways can result in an increase in their stability domain [START_REF] Chauvet | Hierarchical functional organization of formal biological systems: a dynamical approach. I. The increase in complexity by self-association increases the domain of stability of a biological system[END_REF]; this property can be viewed as non-trivial, because an increase in the complexity of artificial systems often results in a decrease of their stability domain.

In our first paper [START_REF] Morillon | Modelling two associated biochemical pathways[END_REF], we showed that the parallel association of two Goodwin metabolic pathways, consisting of a set of nonlinear coupled first-order ordinary differential equations, can increase the stability domain of the unique steady state. Moreover, analytical and numerical studies of the stability of this system show that the association of two units can lead to an increase of the domain of stability. From a biological point of view, the results suggest that the exchange of matter between compartments (e.g. cells, organelles,. . . ) may be a source of stability for the cell metabolism.

The question arises whether this unusual property can be showed in case of nonlinear oscillations. How coupling with exterior pools of metabolites modifies the behaviour of the oscillations of a Goodwin metabolic system.

Moreover the role of spatial heterogeneity for chemical reactions has been the subject of much research (see, for instance, [START_REF] Marmillot | Patterns of spatiotemporal organization in an "ambiquitous" enzyme model[END_REF][START_REF] Capasso | Asymptotic behaviour of reaction-diffusion systems in population and epidemic models. The role of cross diffusion[END_REF] and references therein). Actually the role of exterior pools with respect to the asymptotic behaviour of such systems has been relatively neglected. On the other hand, the presence of exterior pools appears to be rather natural [START_REF] Capasso | Asymptotic behaviour of reaction-diffusion systems in population and epidemic models. The role of cross diffusion[END_REF].

In Section 2 the biological model is presented. In Section 3 some definitions and concepts of the harmonic balancing method are given. In Section 4 we apply this technique in order to compare the oscillations of the systems with or without coupling. In Section 5 numerical experiments are given. In last section, some conclusions are submitted, in particular on stability of steady state. We here investigate the oscillations of the Goodwin metabolic pathway [START_REF] Goodwin | Analytical Physiology of Cells and Developing Organisms[END_REF] by means of the harmonic balancing technique [START_REF] Rapp | Analysis of biochemical phase shift oscillators by a harmonic balancing technique[END_REF], with the addition of some passive diffusion coupling (see Fig. 1). More precisely, we assume that (i) a cell or an organelle u contains a Goodwin metabolic pathway [START_REF] Goodwin | Analytical Physiology of Cells and Developing Organisms[END_REF]; (ii) the metabolite pool P i , i = 1, . . . , n, can exchange matter with the outside pool P * i via passive diffusion, with a non-negative coefficient β i . Let us consider one biological unit which contains a Goodwin metabolic pathway, and apply the laws of mass-balance and enzyme kinetics to this system with passive diffusion between interior and exterior compartments.

P * 1 P * 2 P * 3 P * 4 • • • P * n-1 P * n β 1 ↑↓ β 2 ↑↓ β 3 ↑↓ β 4 ↑↓ β n-1 ↑↓ β n ↑↓ S 0 -→ P 1 -→ P 2 -→ P 3 -→ P 4 -→ • • • -→ P n-1 -→ P n → µ ↑ |
Implementing the transformation of Walter [START_REF] Walter | Kinetic and thermodynamics aspects of biological and biochemical control mechanisms[END_REF][START_REF] Rapp | Analysis of biochemical phase shift oscillators by a harmonic balancing technique[END_REF], the system can be modelled by the following set of differential equations that describe the time evolution of metabolic concentrations :

       dP 1 dt = f(P n ) -b 1 P 1 -β 1 (P 1 -P * 1 ), dP k dt = P k-1 -b k P k -β k (P k -P * k ), k = 2, . . . , n, dP * k dt = β k (P k -P * k ), k = 1, . . . , n, (1)
where P i and P * i are the respective time-dependent concentrations of a given metabolite in the unit and the exterior pools. The coefficient b i caracterizes the positive kinetic constant of the reaction P i → P i+1 in the unit. The reaction function f is given; in the Goodwin metabolic pathway, f describes the allosteric feedback inhibition of the first reaction by the last product, and can be written as:

f(P n ) = 1 1 + (P n ) µ .
where µ is a positive integer. Finally, we assume that association with exterior pools results in passive diffusion, with constant non-negative coefficients β i . If all the β i are zero, the unit is said isolated or non-associated.

Existence and uniqueness of the positive steady state are proved in [START_REF] Morillon | Modelling two associated biochemical pathways[END_REF]. Results on the linear stability of the system at the neighborhood of the unique steady state can be derived by both analytical1 and numerical2 studies for some values of the parameters of the model (see [START_REF] Morillon | Modelling two associated biochemical pathways[END_REF]). Apart the fact that a systematic study of the local properties derived by linearization, may rapidly become cumbersome, the question of discerning whether a system possesses a stable positive steady state or a limit cycle remains to be answered. However, existence of stable oscillations for a Goodwin system without coupling has been proved in [START_REF] Hastings | Existence of periodic solutions for negative feedback cellular control systems[END_REF].

✲ f(P n ) ✲ ✲ G(p) ✛ P n

First order harmonic balancing technique

The harmonic balancing technique for investigating biochemical oscillators can be successfully employed, in particular when the function f of the feedback loop is positive, and monotonically decreasing with his argument.

Definitions

This mathematical method [START_REF] Mees | The describing function matrix[END_REF][START_REF] Gille | Systèmes asservis non linéaires[END_REF][START_REF] Khalil | Nonlinear Systems[END_REF] approximates the output P n and the nonlinear element f(P n ) by a partial sum of their Fourier series :

P n (t) = Re p k=0 x k exp(ikωt) , f(P n (t)) = Re p k=0 y k exp(ikωt) ,
where the coefficients y k verify :

y 0 = 1 2π 2π 0 f Re p k=0 x k exp(ikt) dt, y k = 1 π 2π 0 f Re p l=0 x l exp(ilt) exp(-ikt) dt, 1 ≤ k ≤ p.
Using the loop's input-output equation :

P n = G f(P n )
where G is the linear filter so that :

G f(P n ) = Re ( G(0) y 0 + G(iω) y 1 exp(iω t) + . . . + G(i p ω) y p exp(i p ω t) ) ,
we obtain p + 1 balance equations :

G(i k ω) y k = x k , 0 ≤ k ≤ p .
Then, an approximation to P n (t) is obtained by solving these equations. Since the Fourier series of the function f converges rapidly and G(p) is an extreme low pass filter, the first order approximation :

P n (t) = x 0 + x 1 cos ωt
should be a very good representation of periodic solution for arbitrary dimension n of this system.

Let us define the equivalent gain in the k-th harmonic as F k = y k /x k . This function F k = y k /x k is known as the describing function of the nonlinearity f. This first order approximation leads us to consider the following balance equations :

E G(0) = y 0 x 0 G(0) = 1 , (2)
F G(iω) = y 1 x 1 G(iω) = 1 , (3)
where

y 0 = 1 2π 2π 0 f(x 0 + x 1 cos t) dt, y 1 = 1 π 2π 0 f(x 0 + x 1 cos t) cos t dt.
For single loop control system, E is frequency independent and the zero balance equation [START_REF] Chauvet | Hierarchical functional organization of formal biological systems: a dynamical approach. I. The increase in complexity by self-association increases the domain of stability of a biological system[END_REF] gives us the mean value x 0 as function of the amplitude x 1 of the oscillation. Then, the first balance equation ( 3) becomes a function of the amplitude x 1 and the frequency ω, and it is possible to separate this x 1 -ω relation :

G(iω) = 1 F (x 1 )
.

The G-and 1/F -curves are studied in the complex plane. If they intersect then the balance conditions are satisfied and a periodic solution exists, otherwise the system cannot oscillate. For f monotonically decreasing, it can be proved that 1/F is a continue negative real function, and the set of amplitudes is a connected compact set [START_REF] Rapp | Analysis of biochemical phase shift oscillators by a harmonic balancing technique[END_REF]. An intersection of G(iω)and 1/F (x 1 )-curves will occur only if :

1 F Max ≤ G(iω c ) ≤ 1 F min
where ω c is the frequency so that the G(iω)-curve first intersects the negative real axis. This G (i ω) -curve is often called the Nyquist contour of the linear element.

Functions R and R *

Let us consider the nonlinear system (1) represented by a feedback connection of a linear time-invariant dynamical system and a nonlinear element f(P n ), as shown in Figure 2. We study the behaviour of this autonomous system, and we are interested in investigating the existence of periodic solutions.

Using the simultaneous loop balance conditions, we can construct numerically a function R of the reaction constants where if R > 1 the system oscillates and if R ≤ 1 there is no periodic solution. Comparing the value of R and the numerical solution by Runge-Kutta method via MATLAB software we find that this technique is very accurate.

In our case (see [START_REF] Rapp | Analysis of biochemical phase shift oscillators by a harmonic balancing technique[END_REF] for technical details), the unique condition for an intersection G(iω c ) ≤ 1/F min leads to define R as functions of the parameters b i and β k :

R(b i , β k ) = |F min | b 1 • • • b n • G(iω c ) G(0) (4) so that R(b i , β k ) > 1 is equivalent to this latter condition.
An algorithm is then possible to determine the existence of limit cycle.

In the sequel, let us denote R ≡ R(b i , 0) when the parameters β k are zero, and R * otherwise. Let us note that the nonlinear block is the same for these cases. When we compare an isolated system and an associated system with identical parameters, except β k , we can prove, in special cases, the inequalities R * < R and ω * c ≤ ω c where the star * superscript represents an associated system. Then we observe an increasing stability in the sense that the amplitude and frequency decrease.

For the analysis of the non-associated system, i.e. when all β k = 0, we introduce the linear filter G :

G(p) = n k=1 (p + b k ) -1
where p denotes the differentiation operator. After one or several couplings, the harmonic balance method leads us to consider the following filter :

G * (p) = n k=1 (p + b k )(p + β k ) (p + b k )(p + β k ) + pβ k G(p).
Proposition.

1. G (i ω) = G * (i ω) if and only if ω = 0 or β k = 0 , ∀ k = 1, . . . , n ;
2. The arguments of G(i ω) and G * (i ω) verify the following formulae :

ArgG * (i ω) = ArgG(i ω) -ψ(ω) , ψ(ω) = n k=1 ψ k (ω) , with ψ k (ω) = Arctan (b k + 2 β k ) ω b k β k -ω 2 -Arctan (b k + β k ) ω b k β k -ω 2 ;
(iii) The modulus of G(i ω) and G * (i ω) verify the following estimation :

|G * (i ω)| ≤ |G(i ω)| ∀ ω > 0 .
Proof. The expressions for G(iω) and G * (iω) are :

G(iω) = n k=1 1 (ω 2 + b 2 k ) 1/2 exp -i n k=1 φ k (ω) G * (iω) = G(iω) n l=1 (ω 2 + b 2 l ) 1/2 (ω 2 + β 2 l ) 1/2 (ω 2 + c 2 l ) 1/2 (ω 2 + d 2 l ) 1/2 exp -i n l=1 ψ l (ω) with φ k (ω) = Arctan ω b k , and
ψ l (ω) = Arctan ω c l + Arctan ω d l -Arctan ω b l -Arctan ω β l , if β l = 0 , ψ l (ω) = φ l (ω) if β l = 0 , c l = (b l + 2β l ) -b 2 l + 4β 2 l 2 > 0 , d l = (b l + 2β l ) + b 2 l + 4β 2 l 2 > 0 .
Let us consider the linear filters G and G * with a unique coupling, i.e. with β l = β > 0 and β j = 0 for all j = l. The expressions of G(i ω) and G * (i ω) are deduced of the above, and the equality G(i ω) = G * (i ω) leads to

ω 2 + b 2 l ω 2 + β 2 = ω 2 + c 2 l ω 2 + d 2 l
, and ψ l = 0 , equivalent to ω = 0 or β = 0 .

Moreover, since we have: ψ l (0) = 0, lim β→0 + ψ l (ω) = 0, and b l β = c l d l when ω = 0 or β = 0, the converse is obvious. The arguments of G(i ω) and G * (i ω) verify: ArgG * (i ω) = ArgG(i ω) -ψ l (i ω) and, by means of Arctan property:

ψ l (ω) = Arctan (b l + 2 β) ω b l β -ω 2 -Arctan (b l + β) ω b l β -ω 2 .
The identities:

c 2 l +d 2 l = b 2 l +4 β 2 +2 b l β, and c l d l = b l β lead to the inequality: ∀ ω ∈ R, |G * (i ω)| ≤ |G(i ω)|.
By successive iterations, we obtain the same results with multiple couplings. Corollary. In the complex plane, the G * (i ω)-curve, for ω > 0, is nearer the origin than the G(i ω)curve. Proof. Let us note that we have:

G * (0) = G(0) ∈ R * + , and ∀ ω > 0, ArgG(i ω) = - n k=1 Arctan ω bk < 0.
If ω is large enough, the Arctan property gives us ψ l (ω) < 0, then ArgG * (i ω) > ArgG(i ω). Moreover, as positive ω increases, the modulus and the argument of G(iω) and G * (iω) decrease monotonically:

lim ω→+∞ G(i ω) = lim ω→+∞ G * (i ω) = 0 , lim ω→+∞ ArgG(i ω) = lim ω→+∞ ArgG * (i ω) = -n π 2 .
The curves spiral to origin. Since the G * (i ω)-curve does not intersect the G(i ω)-curve except at ω = 0 or if β k ≡ 0, for all k, and since the modulus and the arguments of G * (i ω) are smaller than those of G(i ω) when ω is large enough, we deduce than the G * (i ω)-curve is nearer the origin in the complex plane than the G(i ω)-curve.

The G(iω)-and G * (iω)-curves first cross the negative real axis at frequencies ω c and ω * c , if they exist, so that :

n i=1 Arctan ω c b i = π (5) n i=1 Arctan ω * c b i + ψ i (ω * c ) = π (6)
where

ψ i (ω * c ) = Arctan (b i + 2β i )ω * c b i β i -(ω * c ) 2 -Arctan (b i + β i )ω * c b i β i -(ω * c ) 2 .
The expression simplifies still further when the following special case is considered. Let us note that a solution ω c > 0 of Eq. 5 exists if and only if n ≥ 3. Proposition.

1. If β l b l > ω 2 c for all l, then ω * c ≤ ω c ; 2. If b k ≡ b and β k ≡ β > 0, 1 ≤ k ≤ n, then (a) ω * c ≤ ω c if and only if β ≥ b tan 2 π n . (b) R * < R . Proof. If β l b l > ω 2
c , for all l = 1, . . . , n, then ψ l (ω c ) > 0 , and

n l=1 Arctan ω c b l + ψ l (ω c ) > π ;
so we obtain:

ω * c ≤ ω c . If b k ≡ b and β k ≡ β > 0, 1 ≤ k ≤ n, we obtain first: ω c = b tan π
n , and we consider the formula 6:

n l=1 Arctan ω * c b + ψ l (ω * c ) = π .
Let us note that ψ l (ω c ) = 0 according to the case β = 0. Let us denote ψ ≡ n ψ 0 with

ψ 0 (ω) = Arctan (b + 2 β) ω b β -ω 2 -Arctan (b + β) ω b β -ω 2 .
We observe that, under this form, ψ 0 (ω) = 0 if and only if ω = 0, or √ b β, or +∞ to the limit. In R + , the function ψ 0 is positive in the interval 0,

√ b β . Let us denote ω * 0 = √ b β. If β ≥ b tan 2 π n , then ∀ ω ∈ 0, √ b β , ψ 0 > 0, and n Arctan ω * 0 b ≥ π, this implies: ω * c ≤ ω * 0 . Since n ψ 0 (ω * c ) ≥ 0, we deduce: n Arctan ω * c b < π, this implies: ω * c ≤ ω c . Conversely, if ω * c ≤ ω c , then ω * c ≤ b tan π n ,
and

n l=1 Arctan ω * c b l = n Arctan ω * c b ≤ π , this leads: ψ (ω * c ) ≥ 0, then ω * c ∈ 0, √ b β . Secondly, it is sufficient to prove |G * (iω * c )| < |G(iω c )| . The representation of G * leads to: |G * (iω)| < |G(iω)| ∀ω > 0 .
ω c = bω * c 5 + b 3 + 2bβ 2 ω * c 3 + 2b 3 β 2 ω * c (b + β)ω * c 4 + (b 2 β + b 3 + bβ 2 )ω * c 2 + b 3 β 2 . Moreover, |G(iω d )| = |G * (iω * c )| implies: ω 2 d = ω * c 4 + (2bβ + 4β 2 )ω * c 2 ω * c 2 + β 2 .
So, we obtain: 2. The concentrations P i , 1 ≤ i ≤ n, of steady states are not modified by coupling; when steady state is unstable, oscillations appear [START_REF] Morillon | Modelling two associated biochemical pathways[END_REF] and coupling decreases the frequency of the oscillations of the system if the diffusion coefficient β is large enough versus kinetic coefficient b. Moreover, if the number of reactants increases, then the stability effect of pools is more effective:

ω 2 d -ω 2 c = C 1 ω * c 12 + C 2 ω * c 10 + C 3 ω * c 8 + C 4 ω * c 6 + C 5 ω * c 4 + C 6 ω * c 2 with C 1 = β 2 + 2 bβ > 0, C 2 = 6 b 3 β + 12 bβ 3 + 6 b 2 β 2 + 4 β 4 > 0, C 3 = (4 b 3 β + 2 bβ 3 )(2 bβ + 4 β 2 ) + 12 b 3 β 3 + 9 b 2 β 4 + 6 b 5 β + 3 b 4 β 2 > 0, C 4 = 2 (b 4 β 2 +b 3 β 3 )(2 bβ+4 β 2 )+2 b 3 β 2 (b 2 β+bβ 2 +b 3 )+10 b 3 β 5 +10 b 5 β 3 +3 b 6 β 2 +2 b 7 β > 0, C 5
n → +∞ =⇒ β ≥ b tan 2 π n → 0 =⇒ ω * c ≤ ω c .
Proposition.

If b k = β k ≡ b, for all k, then 1. n = 4 =⇒ ω * c = ω c ; 2. n > 4 =⇒ ω * c ≤ ω c . Proof. If b k = β k ≡ b , ∀ k, then we have: ψ 0 (ω) = Arctan 2 ω 3 - √ 5 b + Arctan 2 ω 3 + √ 5 b -2 Arctan ω b .
The expression of ψ 0 (ω c ) gives us:

• if n = 4, Arctan 2 3+ √ 5 + Arctan 2 3- √ 5 = π 2 , then ω * c = ω c ; • if n ≥ 4, then ψ 0 (ω c ) = Arctan 3 tanπ/n 1-tan 2 π/n
-2 π n ≥ 0, and the result is obtained:

ω * c ≤ ω c . Proposition.
Let us assume that the systems admit both oscillations, i.e. the G(iω)-and 1/F (x)-curves intersect in the complex plane at point -x of negative real axis (see [START_REF] Rapp | Analysis of biochemical phase shift oscillators by a harmonic balancing technique[END_REF]), and the G * (iω)and 1/F (x)-curves intersect at point -x * < 0 (so, we consider the amplitudes x, x * > 0).

Then: x * < x . Proof. The above results state that, since |G * (iω)| ≤ |G(iω)|, the G * (iω)-curve is nearer the origin in the complex plane than the G(iω)-curve. The G * (i ω) -curve spirals into the origin in a clockwise direction faster than the G (i ω) -curve. This is sufficient to prove the above result. Indeed, when we consider the associated system, the balance method gives us the following first harmonic equation : G * (iω) = 1/F (x) where F is a real function of x, the amplitude of the oscillation, and ω, the frequency (coupling does not modify this function F ). If the G(iω)-, respectively G * (iω)-, and 1/F (x)-curves intersect in the complex plane, oscillations occur with amplitude x, respectively x * , and the amplitude x * of a periodic solution is smaller when coupling exists (the point -x * of the negative real axis is nearer the origin than the point -x).

Numerical experiments

Let us compute the expressions R and R * . First, we need the expressions of G and G * , then the approximate solutions ω c and ω * c of Eqs 5 and 6. Simulations are obtained using the original system:

           dS 1 dτ = α 0 1 + K S µ n -α 1 S 1 -g 1 (S 1 -S * 1 ), dS k dτ = α k-1 S k-1 -α k S k -g k (S k -S * k ), k = 2, . . . , n, dS * k dτ = g k (S k -S * k ), k = 1, . . . , n, (7)
Defining new variables t, P 1 , . . . , P n , P * 1 , . . . , P * n as:

t = α 1 • • • α n K 1/µ α 0 1/n τ ,
P 1 = α -1 0 α 1 • • • α n-1 K 1/µ α 0 1/n S 1 , P k = (α 1 • • • α n-1 α 0 ) -1 α 1 • • • α n-1 K 1/µ α 0 k/n S k , k = 2, . . . , n , P * k = (α 1 • • • α n-1 α 0 ) -1 α 1 • • • α n-1 K 1/µ α 0 k/n S * k .
The resulting equation is equation [START_REF] Capasso | Asymptotic behaviour of reaction-diffusion systems in population and epidemic models. The role of cross diffusion[END_REF] where :

b k = α k α 1 • • • α n-1 K 1/µ α 0 -1/n , β k = g k α 1 • • • α n-1 K 1/µ α 0 -1/n .

Non-associated system

Simulations are obtained, and results are shown in Table 1 using the original system [START_REF] Khalil | Nonlinear Systems[END_REF]. Direct simultations using MATLAB software confirm this fact.

Associated system

Some results are shown in Table 2 where the existence of oscillations are pointed out.

Comparing the value of R * and the numerical solution of the set of equations ( 7), we observe that the values of R * confirm the existence or not of a periodic solution.

Comparison between isolated and associated systems

The results in Table 3 exhibit the effects of coupling on the value R * versus R. When g k are not zero, the numerical values of R * are always smaller than those of R. The stability of associated system is increased.

Frequency comparison

We suppose: α j = a, j = 1, . . . , n, and g k = g, k = 1, . . . , n. Let us note that the condition β > b tan 2 π n for the nondimensioned system (1) is equivalent to the condition g > a tan 2 π n for the original system [START_REF] Khalil | Nonlinear Systems[END_REF]. We can verify the property ω * c < ω c when this condition is satisfied (see, for instance, Table 4 where the opposite cases are shown). In these cases, we have also obtained R * < R. Extensive simulations have conducted us to the special case: if b k ≡ b, β k ≡ β, for all k = 1, . . . , n, so that b = β tan 2 π n , then n = 22 has given ω * c = ω c .

Discussion

Recently, much attention has been paid to the inhomogeneity and anisotropy of the cell milieu. Although the existence of cellular microenvironments has been demonstrated experimentally in many cases, the physiological significance of cellular organization remains open for speculation. Especially, little is known on the consequences of pools existence on the dynamic behaviour of metabolic pathways, while dynamic properties of metabolic pathways in homogeneous media have been studied for several decades (see, for instance, [START_REF] Rapp | Analysis of biochemical phase shift oscillators by a harmonic balancing technique[END_REF]). Since the association of metabolic chains can result in an increase in their stability domain (review in [START_REF] Chauvet | Hierarchical functional organization of formal biological systems: a dynamical approach. I. The increase in complexity by self-association increases the domain of stability of a biological system[END_REF]), we have studied, in the present work, some effects of pools existence on the oscillations of a Yates-Pardee (or Goodwin) metabolic chain. These effects confirm, in part, Chauvet's hypothesis of increasing stability via association of metabolic units.

Our results show that coupling with exterior pools of metabolites modifies the behaviour of the oscillations of a Goodwin metabolic system in the following way : (i) periodic solutions of the system with coupling have lower amplitudes and frequencies; (ii) coupling can give rise to a steady state, instead of a periodic solution. These observations are consistent with the increase of stability observed when two Goodwin-Yates-Pardee metabolic pathways are associated [START_REF] Morillon | Modelling two associated biochemical pathways[END_REF].Thus the stability of the Goodwin-Yates-Pardee metabolic pathway, which is the most commonly encountered form of regulated biochemical pathway, is enhanced when matter can be exchanged with the outside.

Figure 1 :

 1 Figure 1: Goodwin biochemical pathway with couplings

Figure 2 :

 2 Figure 2: Single feedback loop with a nonlinear function f and a linear filter G(p)

  The functions |G * (iω)| and |G(iω)| versus ω are continuous and monotonically decreasing with increasing positive frequency.Let us denote ω d so that|G * (iω * c )| = |G(iω d )| . If ω c < ω d , |G(iω d )| < |G(iω c )|, then |G * (iω * c )| < |G(iω c )| . It remains to prove: ω c < ω d , so |G * (iω * c )| < |G(iω c )| is always verified. The equality:

= 5 b 6 β 4 + 10 b 5 β 5 + 4 b 7 β 3 >Applications 1 .

 4531 0 and C 6 = 2 b 7 β 5 . Lastly, we deduce: ω c < ω d , then R * < R . According to the condition: β ≥ b tan 2 π n , when n ≥ 3, it is possible to find ω c and ω * c solutions of Eqs 5 and 6. This condition leads to β ≥ +∞ when n = 2, β ≥ 3 b when n = 3, β ≥ b/3 when n = 6 . In the case n = 2, when b k ≡ b and β l ≡ β, the formulae 5 and 6 have no solutions: no oscillations occur. In the cases n > 2, couplings with large β give smaller frequencies.

Table 1 :

 1 Nonassociated system: α j = 1, j = 1, . . . , n -1, µ = 5

	n	α 0	K	α n	R	Oscillation
	5	50	0.01	0.5	0.8673	No
	5	50	1	1	0.9851	No
	7	5	0.1	0.5	0.9730	No
	7	50	1	5	1.1115	Yes
	8	5	0.1	0.5	1.0811	Yes
	8	5	0.1	1	0.8653	No

Table 2 :

 2 Associated system: α j = 1, j = 1, . . . , n -1, g l = 1, l = 1, . . . , n, µ = 3

	n	α 0	K	α n	R *	Oscillation
	7	50	0.1	0.5	0.8735	No
	8	10	10	1	0.9953	No
	9	150	0.1	0.5	1.2199	Yes
	10	50	1	0.5	1.3253	Yes
	12	1	2	0.5	0.8911	No
	15	200	0.01	0.1	1.2414	Yes

Table 3 :

 3 Isolated and associated systems: α j = 1, j = 1, . . . , n -1, g l = 1, l = 1, . . . , n, µ = 3

	n	α 0	K	α n	R	R *
	7	50	0.1	0.5	1.2949	0.8735
	8	5	0.1	1	0.8653	0.5992
	9	150	1	1	1.6742	1.2231
	10	10	2	0.1	0.9513	0.7967
	12	100	0.01	5	1.3426	1.0568

Table 4 :

 4 Frequency comparison: α j = a, j = 1, . . . , n, g l = g, l = 1, . . . , n, K = 2

	n	α 0	µ	g	a tan 2 π n	ω c	ω * c
	3	50	6	0.065	0.783	0.452	0.556
	4	200	7	0.046	0.461	0.461	0.502
	4	50	5	0.726	0.363	0.363	0.222
	6	5	6	0.210	0.280	0.486	0.516
	6	200	2	1.020	0.170	0.294	0.155
	8	100	4	1.009	0.086	0.209	0.107

Using the criterion of Liénard and Chipart.

By computing the eigenvalues of jacobian matrix with scientific software such as MATLAB.


Charlotte Brennan
email: charlotte.brennan@wits.ac.za

Arnold Knopfmacher
email: arnold.knopfmacher@wits.ac.za

Descent variation of samples of geometric random variables

Keywords: words, q-enumeration, generating functions, up-down words, Carlitz words, descents

In this paper, we consider random words ω1ω2ω3 • • • ωn of length n, where the letters ωi ∈ N are independently generated with a geometric probability such that P{ωi = k} = pq k-1 where p + q = 1. We have a descent at position i whenever ωi+1 < ωi. The size of such a descent is ωi -ωi+1 and the descent variation is the sum of all the descent sizes for that word. We study various types of random words over the infinite alphabet N, where the letters have geometric probabilities, and find the probability generating functions for descent variation of such words.

Introduction

Recently Cristea and Prodinger [START_REF] Cristea | q-enumeration of words by their total variation[END_REF] considered words ω = ω 1 ω 2 . . . ω n over the infinite alphabet N, where the letters occur independently with the geometric probability P(ω i = j) = pq j-1 , for all n, i, j ∈ N with i ≤ n, where 0 < p, q < 1 and p + q = 1. They defined the total variation to be

u(ω) = n-1 i=1 |ω i+1 -ω i |.
They investigated the probability generating function of a word of given length n with a total variation u(ω) = r.

Three different families of words over N were studied:

• the general case, namely all words where the letters occur with a geometric probability with no restriction,

• words with an up-down pattern <≥<≥ • • • <≥,

• Carlitz words where any two adjacent letters must be distinct.

Previously, several authors had studied the statistic total variation in other contexts. For more information, see the paper [START_REF] Mansour | Enumeration of words by the sum of differences between adjacent letters[END_REF] by Mansour. In our paper, we also consider words over N where the letters are independently generated with the same geometric probability. However, we study the statistic "descent variation". We say we have a descent at position i whenever ω i+1 < ω i . The size of such a descent is the difference ω i -ω i+1 . We study the sum of all such descent sizes for the whole word. We call this sum, the descent variation. We concentrate our study on the same three families of words as described above. In Section 2, we look at the general case. In Section 3, we look at the word with an up-down pattern <≥<≥ • • • <≥ by first considering word of odd length then words of even length. Finally, in Section 4, we look at the descent variation of Carlitz words, where no two adjacent letters have equal value.

We use the "adding a slice" method to find the generating functions, this method has been used in [START_REF] Knopfmacher | On Carlitz compositions[END_REF]. The methodology is that we start with a word of length n ending with a letter with a specified value i and add an extra letter of value j at the end of the word. This leads to functional equations which can be solved to deduce the generating functions.

2 The general case: geometrically distributed words over N Let W n denote the set of words of length n with alphabet N and where n ≥ 1.

For ω ∈ W n , n ≥ 1, ω = ω 1 ω 2 . . . ω n , we define i) a descent at position i whenever ω i > ω i+1 ,
ii) a descent size, as the difference ω i -ω i+1 whenever there is a descent and iii) the descent variation of the word ω, as the sum of all the descent sizes.

Our aim is to find, for given n ≥ 1 and r ≥ 0, the probability that a word ω ∈ W n has descent variation u(ω) = r.

We introduce the following generating functions. Let

f (z, x) = n≥1,r≥0 p(n, r)z n x r , (2.1)
where p(n, r) is the probability that a word of length n has descent variation r. Moreover, let

F (z, x, t) = n≥1,r≥0,i≥1 g(n, r, i)z n x r t i , (2.2)
where g(n, r, i) is the occurrence probability of a word of length n with last letter i (ω n = i) and descent variation r. We have f (z, x) = F (z, x, 1).

In order to find the generating function of words with n letters and descent variation r, we proceed as follows.

Suppose a random word ω of length n has an arbitrary last letter with fixed value i and descent variation r, ω = ω 1 ω 2 . . . ω n-1 i, where r ≥ 0 and i ≥ 1. We add a letter, j, at the end of the word ω, after the last letter i. A new word ω ′ of length n + 1 is obtained,

ω ′ = ω 1 ω 2 . . . ω n-1 ij.
We have two cases depending on the value of j relative to i:

• if j > i, there is no new descent and therefore the descent variation remains the same r, or

• if j ≤ i, we have one more descent, thus the descent variation is now r + i -j.

The following two cases with the expression (2.3) represent the changes that adding the extra letter induces in the generating function F (z, x, t).

i-1 j=1 pq j-1 x i-j t j + j≥i pq j-1 t j = pt x -qt x i + pt (x -1) (x -qt) (1 -qt) (qt) i . (2.3)
According to the "adding a slice" method, see [START_REF] Knopfmacher | On Carlitz compositions[END_REF], the process of going from a word of length n with last letter i to a word of length n + 1 with last letter j is therefore achieved by making the following substitution in our generating function

t i -→ pt x -qt x i + pt (x -1) (x -qt) (1 -qt) (qt) i . (2.4)
Let a n (x, t), n ≥ 1, be the function with the property

F (z, x, t) = n≥1 a n (x, t)z n . (2.5)
For given n ≥ 1, using the substitution given in (2.4) we obtain

a n+1 (x, t) = pt x -qt a n (x, x) + pt (x -1) (x -qt) (1 -qt) a n (x, qt), (2.6)
which, after multiplication by z n+1 and then summing up for n ≥ 1, leads to

F (z, x, t) = a 1 (x, t)z + z pt x -qt F (z, x, x) + z pt (x -1) (x -qt) (1 -qt) F (z, x, qt), (2.7)
where a 1 (x, t) = k≥1 pq k-1 t k = pt 1-qt . This can now be solved using the following elementary lemma:

Lemma 1 A recursion of the form F (z, t) = f (z, t) + g(z, t)F (z, tq) has solution F (z, t) = j≥1 f (z, tq j-1 ) j-1 i=1 g(z, tq i-1
).

If we use Lemma 1 with

f (z, t) = ptz 1 -qt + ptz x -qt F (z, x, x) and g(z, t) = ptz (x -1) (x -qt) (1 -qt) ,
we obtain

F (z, x, t) = j≥1 ptzq j-1 1 -tq j + ptzq j-1 x -tq j F (z, x, x) • j-1 i=1 ptzq i-1 (x -1) (x -tq i )(1 -tq i )
.

(2.8)

For t = x, we get F (z, x, x) = j≥1 pxzq j-1 1 -xq j + pzq j-1 1 -q j F (z, x, x) • j-1 i=1 pzq i-1 (x -1) (1 -q i )(1 -xq i ) ,
and thus

F (z, x, x) = j≥1 (zp) j x(x-1) j-1 q ( j 2 )
(q;q)j-1(xq;q)j 1 -j≥1 (zp) j (x-1) j-1 q ( j 2 )

(q;q)j (xq;q)j-1

, where (x; q) n := n-1 j=0 (1 -xq j ) for n ≥ 1, and (x; q) 0 = 1. Substituting the above expression for F (z, x, x) into (2.8) we obtain

F (z, x, t) = j≥1 pztq j-1 1 -tq j + zptq j-1 x -tq j k≥1 (zp) k x(x-1) k-1 q ( k 2 ) (q;q) k-1 (xq;q) k 1 -k≥1 (zp) k (x-1) k-1 q ( k 2 ) (q;q) k (xq;q) k-1 • zpt(x -1) j-1 q ( j-1 2 )
x j-1 ( tq x ; q) j-1 (tq; q) j-1

.

For t = 1 we obtain

f (z, x) = F (z, x, 1) = j≥1 pzq j-1 1 -q j + zpq j-1 x -q j k≥1 (zp) k x(x-1) k-1 q ( k 2 ) (q;q) k-1 (xq;q) k 1 -k≥1 (zp) k (x-1) k-1 q ( k 2 ) (q;q) k (xq;q) k-1 • zp(x -1) j-1 q ( j-1 2 )
x j-1 ( q x ; q) j-1 (q; q) j-1

.

(2.9)

We can test our formula, by putting x = 1:

f (z, 1) = F (z, 1, 1) = z 1 -z ,
which corresponds to the case when neither the last letter nor the descent variation of the words are considered. Furthermore,

∂f (x, z) ∂x x=1 = qz 2 (1 -q 2 )(1 -z) 2 .
Thus the mean descent variation of a word with n letters is

[z n ] ∂f (x, z) ∂x x=1 = [z n ] qz 2 (1 -q 2 )(1 -z) 2 = q 1 -q 2 (n -1).
Cristea and Prodinger in [START_REF] Cristea | q-enumeration of words by their total variation[END_REF] found twice that value for the mean total variation, which is to be expected by symmetry. However, the variance of the descent variation and of the total variation are quite different. We obtain

∂ 2 f (x, z) ∂x 2 x=1 = 2q 2 z 2 -1 + q(-2 + z) -q 2 (2 -2z + z 2 ) + q 3 (-1 + z) + z -z 2 (1 -q 2 ) 2 (1 + q + q 2 )(-1 + z) 3 ,
where

[z n ] ∂ 2 f (x, z) ∂x 2 x=1 =      q 2 n 2 (1+q+q 2 )+n(-3+q-q 2 +2q 3 )-2(-2+q-q 2 +q 3 ) (1-q 2 ) 2 (1+q+q 2 )
, for n > 1 0, for n = 1.

The variance for the descent variation is therefore

nq q 4 + 3q 2 + 1 (q 2 -1) 2 (q 2 + q + 1) - q q 2 -q + 1 2 (q 2 -1) 2 (q 2 + q + 1) , for n > 1.
3 Up-down words over N following the pattern <≥<≥ • • • <≥

Up-down words of odd length

The aim of this subsection is to find the generating function of up-down words with pattern <≥<≥ • • • <≥ with odd length and given descent variation. Corresponding generating functions for such words by total variation and by the number of rises have been found in [START_REF] Cristea | q-enumeration of words by their total variation[END_REF][START_REF] Cristea | q-enumeration of up-down words by number of rises[END_REF]. Note, these are not probability generating functions as we still need to normalize the coefficients of z 2n+1 by dividing by the total probability of having alternating words of length 2n + 1.

The methodology is the same as for the general case in the previous section, we introduce the following generating functions. Let

f (z, x) = n≥0,r≥0 p(2n + 1, r)z 2n+1 x r , (3.1)
where p(2n + 1, r) is the probability that a word of length 2n + 1 has descent variation r, and let

F (z, x, t) = n≥0,r≥0,i≥1 g(2n + 1, r, i)z 2n+1 x r t i = n≥0 a 2n+1 (x, t)z 2n+1 (3.2)
where g(2n+1, r, i) is the occurrence probability of a word of length 2n+1, with last letter i (ω 2n+1 = i) and descent variation r. Again, f (z, x) = F (z, x, 1). We use the "adding a slice" method but this time we need to add two letters at a time. We proceed inductively. Assume that we have a random word ω = ω 1 . . . ω 2n i, of odd length 2n + 1 with last letter i, i.e., w 2n+1 = i and w 2n ≥ i and descent variation r, with arbitrarily fixed i ≥ 1 and r > 0. We add two letters, k and j, after the last letter i. A new word, w ′ of odd length 2n + 3 is obtained,

ω ′ = ω 1 . . . ω 2n ikj.
To carry on the up-down pattern, we need i < k, and k ≥ j. The descent variation of ω ′ is now r + k -j.

The following expression summarises the situation

k>i pq k-1 j≤k pq j-1 • t j • x k-j = p 2 xt(xq) i (x -qt)(1 -qx) - p 2 t 2 q(q 2 t) i (x -qt)(1 -q 2 t) . (3.3)
Again, using the "adding a slice" method, this leads to the following substitution

t i -→ p 2 t x -qt • x 1 -qx (qx) i - qt 1 -q 2 t (q 2 t) i , (3.4)
or, equivalently,

t i -→ A(x, t) x (1 -qx) (qx) i -B(t) • (q 2 t) i , (3.5)
where A(x, t) = p 2 t x-qt and B(t) = qt 1-q 2 t . We therefore obtain, for n ≥ 0,

a 2n+3 (x, t) = A(x, t) x 1 -qx a 2n+1 (x, qx) -B(t) a 2n+1 (x, q 2 t) , (3.6)
which, after multiplication by z 2n+3 and then summing up for n ≥ 0, leads to

F (z, x, t) = a 1 (x, t)z + A(x, t)z 2 x 1 -qx F (z, x, qx) -B(t)F (z, x, q 2 t) , (3.7)
with a 1 (x, t) = k≥1 pq k-1 t k = pt 1-qt . We solve the functional equation using Lemma 1, where q is replaced by q 2 . Thus

F (z, x, t) = k≥1 ptq 2k-2 z 1 -tq 2k-1 + A(x, tq 2k-2 ) z 2 x 1 -qx F (z, x, qx) (-1) k-1 z 2(k-1) k-1 i=1 A(x, tq 2i-2 ) B(tq 2i-2 ).
In order to obtain F (z, x, qx), let t = qx F (z, x, qx)

= k≥1 pxq 2k-1 z 1 -xq 2k + A(x, xq 2k-1 ) z 2 x 1 -qx F (z, x, xq) (-1) k-1 z 2(k-1) k-1 i=1 A(x, xq 2i-1 )B(xq 2i-1 )
which yields

F (z, x, qx) = k≥1 pxq 2k-1 z 2k-1 (-1) k-1 1 -xq 2k k-1 i=1 A(x, xq 2i-1 )B(xq 2i-1 ) 1 - k≥1 A(x, xq 2k-1 )x(-1) k-1 z 2k 1 -qx k-1 i=1 A(x, xq 2i-1 )B(xq 2i-1
) .

We have

k-1 i=1 A(x, xq 2i-1 ) = k-1 i=1 p 2 q 2i-1 1 -q 2i = p 2(k-1) q (k-1) 2 (q 2 ; q 2 ) k-1 and k-1 i=1 B(xq 2i-1 ) = k-1 i=1 xq 2i 1 -xq 2i+1 = x k-1 q k(k-1) (xq 3 ; q 2 ) k-1 ,

and thus

k-1 i=1

A(x, xq 2i-1 )B(xq 2i-1 ) = p 2(k-1) q (2k-1)(k-1) x k-1 (q 2 ; q 2 ) k-1 (xq 3 ; q 2 ) k-1 .

Moreover, since

k-1 i=1 A(x, tq 2(i-1) ) = k-1 i=1 p 2 tq 2(i-1) (x -tq 2i-1 ) = (p 2 t) k-1 q (k-2)(k-1) x k-1 ( tq x ; q 2 ) k-1 and k-1 i=1 B(tq 2(i-1) ) = k-1 i=1 tq 2i-1 1 -tq 2i = t k-1 q (k-1) 2 (tq 2 ; q 2 ) k-1 , we have k-1 i=1 A(x, tq 2(i-1) )B(tq 2(i-1) ) = (p 2 t 2 ) k-1 q (2k-3)(k-1) x k-1 ( tq x ; q 2 ) k-1 (tq 2 ; q 2 ) k-1
.

These finally yield an expression for F (z, x, t):

F (z, x, t) = n≥1 p 2 tq 2n-2 z 2 x (x -tq 2n-1 )(1 -qx) • k≥1 (-1) k-1 p 2k-1 q k(2k-1) x k z 2k-1 (1-xq 2k )(q 2 ;q 2 ) k-1 (xq 3 ;q 2 ) k-1 1 -k≥1 (-1) k-1 p 2k q k(2k-1) x k z 2k (q 2 ;q 2 ) k (qx;q 2 ) k + ptq 2(n-1) z 1 -tq 2n-1 × (-1) n-1 (p 2 t 2 z 2 ) n-1 q (2n-3)(n-1)
x n-1 ( qt x ; q 2 ) n-1 (tq 2 ; q 2 ) n-1

.

Setting t = 1 gives the sought generating function for the descent variation,

f (z, x) = F (z, x, 1) = n≥1 p 2 q 2n-2 z 2 x (x -q 2n-1 )(1 -qx) • k≥1 (-1) k-1 p 2k-1 q k(2k-1) x k z 2k-1 (1-xq 2k )(q 2 ;q 2 ) k-1 (xq 3 ;q 2 ) k-1 1 -k≥1 (-1) k-1 p 2k q k(2k-1) x k z 2k (q 2 ;q 2 ) k (qx;q 2 ) k + pq 2(n-1) z 1 -q 2n-1 × (-1) n-1 (p 2 z 2 ) n-1 q (2n-3)(n-1)
x n-1 ( q x ; q 2 ) n-1 (q 2 ; q 2 ) n-1

.

(3.8)

In order to compute probabilities, the coefficients of z 2n+1 in f (z, x) must be normalized by dividing by the coefficients of z 2n+1 in f (z, 1). Remarks As a simple check of (3.8), a series expansion of f (z, 1) and f (z, x), for p = q = 1/2 give 2/7 as the probability that a three letter word is up-down, and 34/63 for the mean descent variation. These values are easily verified by direct calculation. For example, the mean descent variation in the case of a three letter word is evaluated by the sum

∞ i=1 pq i-1 ∞ k=i+1 pq k-1 k j=1 (k -j)pq j-1 .
By contrast, direct calculation gives a value of 8/7 for the mean total variation in this case, for p = q = 1/2. However, this does not correspond to a series expansion of the generating function in [START_REF] Cristea | q-enumeration of words by their total variation[END_REF], indicating that the result there contains errors. Similarly, we find that the total variation generating function for up-down words of even length in [START_REF] Cristea | q-enumeration of words by their total variation[END_REF] is incorrect.

Up-down words of even length

In this subsection, we are finding the generating functions of up-down words of even length with a given descent variation.

Assume we have a random word ω ′ = ω 1 . . . ω 2n+1 l of even length 2n + 2, with last letter l, and descent variation r, with arbitrarily fixed l ≥ 1 and r > 0. We take a word of odd length 2n + 1 ending with j and add the letter l. We use and adapt the result we had for a word of odd length in (3.3). Then according to the pattern, we need l > j, and therefore the descent variation of ω ′ is r + l -j. Thus we have the following expression k>i p q k-1 j≤k p q j-1 x k-j l>j p q l-1 z t l .

(3.9)

We introduce the generating function

g(z, x) := n≥0,r≥0 p(2n + 2, r)z 2n+2 x r , (3.10)
where p(2n + 2, r) is the probability that an up-down word of length 2n + 2 has descent variation r.

Since l>j pq l-1 zt l = pzt 1-qt (qt) j , we can obtain the sum (3.9) by replacing t j in the sum (3.3) by pzt 1-qt (qt) j . Therefore, for t = 1 we simply use F (z, x, t) in the case of up-down words of odd length, replace t by q and multiply by z to get g(z, x) = zF (z, x, q), i.e., the generating function for descent variation in up-down words of even length is

g(z, x) = n≥1 p 2 q 2n-1 z 3 x (x -q 2n )(1 -qx) • k≥1 (-1) k-1 p 2k-1 q k(2k-1) x k z 2k-1 (1-xq 2k )(q 2 ;q 2 ) k-1 (xq 3 ;q 2 ) k-1 1 -k≥1 (-1) k-1 p 2k q k(2k-1)
x k z 2k (q 2 ;q 2 ) k (qx;q 2 ) k

+ pq 2n-1 z 2 1 -q 2n × (-1) n-1 (p 2 q 2 z 2 ) n-1 q (2n-3)(n-1)
x n-1 ( q 2 x ; q 2 ) n-1 (q 3 ; q 2 ) n-1

.

Once again, to compute the probabilities we need to divide [z 2n ]f (z, x) by [z 2n ]f (z, 1).

Carlitz words

In this section, we consider a special type of words called Carlitz words. These are are words over the alphabet N, where the letters again have geometric probabilities and with the property that any two adjacent letters are distinct. In this section, we find the generating function for a Carlitz word with given descent variation. Various people have studied the related topic of Carlitz compositions of natural numbers, namely, Carlitz [START_REF] Carlitz | Restricted compositions[END_REF], Knopfmacher and Prodinger [START_REF] Knopfmacher | On Carlitz compositions[END_REF], Goh and Hitczenko [START_REF] Goh | Average Number of Distinct Path Sizes in a Random Carlitz Composition[END_REF], and Louchard and Prodinger [START_REF] Louchard | Probabilistic Analysis of Carlitz Compositions[END_REF].

We use the same notation as for the general case, in Section 2, namely

f (z, x) = n≥1,r≥0 p(n, r)z n x r , (4.1)
where p(n, r) is in this case the probability that a Carlitz word of length n has descent variation r, and

F (z, x, t) = n≥1,r≥0,i≥1 g(n, r, i)z n x r t i , (4.2)
where g(n, r, i) is the probability that a Carlitz word of length n, with last letter i (ω n = i) has descent variation r. We have f (z, x) = F (z, x, 1). Again, to obtain probability generating functions we need to normalize the coefficients of z n by dividing by the probability of having a Carlitz word of length n.

We perform the "adding a slice" technique again. Suppose we have a random Carlitz word of length n with last letter i, ω = ω 1 . . . ω n-1 i, with descent variation r, with arbitrarily fixed i ≥ 1 and r ≥ 0. We add a letter, j, after the last letter i, and thus get a new word of length n + 1, ω ′ = ω 1 . . . ω n-1 ij.

The letter j can be either greater than i or less than i, but not equal to i.

• If j > i then the descent variation remains the same.

• If j < i then the new word ω ′ has descent variation r + i -j.

This is reflected by the following expression

i-1 j=1 p q j-1 x i-j t j + j≥i+1 p q j-1 t j . (4.3) Thus, we have the substitution

t i -→ pt (x -qt) x i + pt 1 -qt - px q(x -qt) (qt) i . (4.4)
Again, we let F (z, x, t) = n≥1 a n (x, t)z n . For n ≥ 1, we get from (4.4),

a n+1 (x, t) = pt (x -qt) a n (x, x) + pt 1 -qt - px q(x -qt) a n (x, qt) (4.5)
which, after multiplication by z n+1 and then summing over n ≥ 1, we have

F (z, x, t) = a 1 (x, t)z + pt (x -qt) zF (z, x, x) + pt 1 -qt - px q(x -qt) zF (z, x, qt), (4.6)
where a 1 (x, t) = k≥1 pq k-1 t k = pt 1-qt . Using Lemma 1, with

f (z, x, t) = ptz 1 -qt + ptz x -qt F (z, x, x), and
g(z, x, t) = ptz 1 -qt - pxz q(x -qt) we obtain F (z, x, t) = j≥1 ptq j-1 z 1 -tq j + ptq j-1 z x -tq j F (z, x, x) (pz) j-1 • j-1 i=1 tq i-1 1 -tq i - x q(x -tq i ) . (4.7)
For t = x, we get

F (z, x, x) = j≥1 pxq j-1 z 1 -xq j + pq j-1 z 1 -q j F (z, x, x) (pz) j-1 • j-1 i=1 xq i-1 1 -xq i - 1 q(1 -q i ) (4.8)
and thus F (z, x, x) = j≥1 xq j-1 (pz) j 1-xq j j-1 i=1

xq i-1

1-xq i -1 q(1-q i ) 1 -j≥1 q j-1 1-q j (pz) j j-1 i=1 xq i-1 1-xq i -1 q(1-q i ) . This yields

F (z, x, t) = n≥1   tq n-1 1 -tq n + tq n-1 x -tq n k≥1 xq k-1 (pz) k 1-xq k k-1 i=1 xq i-1 1-xq i - 1 q(1-q i ) 1 -k≥1 q k-1 1-q k (pz) k k-1 i=1 xq i-1 1-xq i - 1 q(1-q i )   × (pz) n • n-1 i=1 tq i-1 1 -tq i - x q(x -tq i )
. (4.9)

After substituting t = 1 we get our desired probability generating function

f (z, x) = F (z, x, 1) = n≥1   q n-1 1 -q n + q n-1 x -q n k≥1 xq k-1 (pz) k 1-xq k k-1 i=1 xq i-1 1-xq i - 1 q(1-q i ) 1 -k≥1 q k-1 1-q k (pz) k k-1 i=1 xq i-1 1-xq i - 1 q(1-q i )   × (pz) n • n-1 i=1 q i-1 1 -q i - x q(x -q i ) . (4.10)
Furthermore, we obtain after substituting x = 1 and t = 1 in (4.7):

f (z, 1) = F (z, 1, 1) = j≥1 pq j-1 z 1 -q j + pq j-1 z 1 -q j F (z, 1, 1) (pz) j-1 • j-1 i=1 q i-1 1 -q i -1 q(1 -q i ) = j≥1 1 1 -q j + 1 1 -q j F (z, 1, 1) (-1) j-1 (pz) j , which yields f (z, 1) + 1 = 1 1 -j≥1 (-1) j-1 (pz) j 1-q j , (4.11) which is the generating function for Carlitz words. (i) In order to compute probabilities, the coefficients of z n in f (z, x) must be normalized by dividing by the coefficients of z n in f (z, 1). (i) The corresponding f (z, 1) for Carlitz words in [START_REF] Cristea | q-enumeration of words by their total variation[END_REF] does not match (4.11), so is incorrect.


L Sunil Chandran

Rogers Mathew
email: rogers@csa.iisc.ernet.in

Balakrishnan

Bipartite Powers of k-chordal Graphs

Keywords: k-chordal graph, hole, chordality, graph power, bipartite power

Let k be an integer and k ≥ 3. A graph G is k-chordal if G does not have an induced cycle of length greater than k. From the definition it is clear that 3-chordal graphs are precisely the class of chordal graphs. Duchet proved that, for every positive integer m, if G m is chordal then so is G m+2 . Brandstädt et al. in [Andreas Brandstädt, Van Bang Le, and Thomas Szymczak. Duchet-type theorems for powers of HHD-free graphs. Discrete Mathematics, 177(1-3): 9-16, 1997.] showed that if G m is k-chordal, then so is G m+2 .

Powering a bipartite graph does not preserve its bipartitedness. In order to preserve the bipartitedness of a bipartite graph while powering Chandran et al. introduced the notion of bipartite powering. This notion was introduced to aid their study of boxicity of chordal bipartite graphs. The m-th bipartite power G [m] of a bipartite graph G is the bipartite graph obtained from G by adding edges (u, v) where dG(u, v) is odd and less than or equal to m. Note that G [m] = G [m+1] for each odd m.

In this paper we show that, given a bipartite graph G, if G is k-chordal then so is G [m] , where k, m are positive integers with k ≥ 4.

Introduction

A hole is a chordless (or an induced) cycle in a graph. The chordality of a graph G, denoted by C(G), is defined to be the size of a largest hole in G, if there exists a cycle in G. If G is acyclic, then its chordality is taken as 0. A graph G is k-chordal if C(G) ≤ k. In other words, a graph is k-chordal if it has no holes with more than k vertices in it. Chordal graphs are exactly the class of 3-chordal graphs and chordal bipartite graphs are bipartite, 4-chordal graphs. k-chordal graphs have been studied in the literature in [START_REF] Bodlaender | Treewidth for graphs with small chordality[END_REF], [START_REF] Sunil Chandran | On the number of minimum cuts in a graph[END_REF], [START_REF] Sunil Chandran | Graphs of low chordality[END_REF], [START_REF] Dourisboure | Compact routing schemes for generalised chordal graphs[END_REF], [START_REF] Feodor | Estimating all pairs shortest paths in restricted graph families: a unified approach[END_REF] and [START_REF] Spinrad | Finding large holes[END_REF]. For example, Chandran and Ram [START_REF] Sunil Chandran | On the number of minimum cuts in a graph[END_REF] proved that the number of minimum cuts in a k-chordal graph is at most (k+1)n Paulraja [START_REF] Balakrishnan | Powers of chordal graphs[END_REF] proved that odd powers of chordal graphs are chordal. Chang and Nemhauser [START_REF] Chang | The k-domination and k-stability problems on sun-free chordal graphs[END_REF] showed that if G and G 2 are chordal then so are all powers of G. Duchet [START_REF] Duchet | Classical perfect graphs[END_REF] proved a stronger result which says that if G m is chordal then so is G m+2 . Brandstädt et al. in [START_REF] Brandstädt | Duchet-type theorems for powers of hhd-free graphs[END_REF] showed that if G m is k-chordal then so is G m+2 , where k ≥ 3 is an integer. Studies on families of graphs that are closed under powering can also be seen in the literature. For instance, it is known that interval graphs, proper interval graphs [START_REF] Raychaudhuri | On powers of interval and unit interval graphs[END_REF], strongly chordal graphs [START_REF] Lubiw | γ-free matrices[END_REF], circular-arc graphs [START_REF] Raychaudhuri | On powers of strongly chordal and circular graphs[END_REF] [START_REF] Flotow | On powers of circular arc graphs and proper circular arc graphs[END_REF], cocomparability graphs [START_REF] Flotow | On powers of m-trapezoid graphs[END_REF] etc. are closed under taking powers.

Subclasses of bipartite graphs, like chordal bipartite graphs, are not closed under powering since the m-th power of a bipartite graph need not be even bipartite. Chandran et al. in [START_REF] Chandran | Chordal bipartite graphs with high boxicity[END_REF] introduced the notion of bipartite powering to retain the bipartitedness of a bipartite graph while taking power. The m-th bipartite power G [m] of a bipartite graph G is the bipartite graph obtained from G by adding edges (u, v) where d G (u, v) is odd and less than or equal to m. Note that G [m] = G [m+1] for each odd m. It was shown in [START_REF] Chandran | Chordal bipartite graphs with high boxicity[END_REF] that the m-th bipartite power of a tree is chordal bipartite. The intention there was to construct chordal bipartite graphs of high boxicity. The fact that the chordal bipartite graph under consideration was obtained as a bipartite power of a tree was crucial for proving that its boxicity was high. Since trees are a subclass of chordal bipartite graphs, a natural question that came up was the following: is it true that the m-th bipartite power of every chordal bipartite graph is chordal bipartite? In this paper we answer this question in the affirmative. In fact, we prove a more general result.

Observation 2. H is an induced subgraph of H

′ . Observation 3. Let k be an integer such that k ≥ 4. If H is k-chordal, then so is H ′ .
Proof: Any hole of size greater than 4 in H ′ cannot have more than one vertex from the same bag, say B v , as such vertices have the same neighborhood. Hence, the vertices of a hole (of size greater than 4) in H ′ belong to different bags and thus there is a corresponding hole of the same size in H.

Theorem 4. Let m, k be positive integers with k ≥ 4. Let G be a bipartite graph. If G is k-chordal, then so is G [m] .
Proof: We prove this by contradiction. Let p denote the size of a largest induced cycle, say C = u 0 u 1 . . . u p-1 u 0 , in G [m] . Assume p > k. Then, p ≥ 6 (since k ≥ 4 and G [m] is bipartite). Between each u i-1 and u i , where i ∈ {0, . . . , p -1}, there exists a shortest path of length not more than m in G (i) . Let P i be one such shortest path between u i-1 and u i in G.

Let H be the subgraph induced on the vertex set

p-1 i=0 V (P i ) in G.
As mentioned in the beginning of this section, construct a graph H ′ from H, where H ′ ∈ B(H), in the following manner: for each

v ∈ V (H), let |B v | = |{P i | 0 ≤ i ≤ p -1, v ∈ V (P i )}| i.e.
, let B v have as many vertices as the number of paths in {P 0 . . . P p-1 } that share vertex v in H. For each i ∈ {0, . . . , p -1}, let Q ′ i = u i-1 Q i be a shortest path between u i-1 and u i in H ′ such that no two paths Q i and Q j (where i = j) share a vertex (i) . From our construction of H ′ from H it is easy to see that such paths exist. Let

Q i = v i,1 v i,2 . . . v i,ri u i , where r i = ||Q i || ≥ 0. Thus, Q ′ i = u i-1 v i,1 v i,2 . . . v i,ri u i . Clearly, ||Q ′ i || = ||P i || ≤ m.
The reader may also note that the cycle C (= u 0 u 1 . . . u p-1 u 0 ) which is present in G [m] will be present in H [m] and thereby in H ′[m] too.

In order to prove the theorem, it is enough to show that there exists an induced cycle of size at least p in H ′ . Then by combining Observation 3 and the fact that H is an induced subgraph of G, we get k ≥ C(G) ≥ C(H) ≥ C(H ′ ) ≥ p contradicting our assumption that p > k. Hence, in the rest of the proof we show that C(H ′ ) ≥ p.

Consider the following drawing of the graph H ′ . Arrange the vertices u 0 , u 1 , . . . , u p-1 in that order on a circle in clockwise order. Between each u i-1 and u i on the circle arrange the vertices v i,1 , v i,2 , . . . , v i,ri in that order in clockwise order. Recall that these vertices are the internal vertices of path Q ′ i . Claim 4.1. In this circular arrangement of vertices of H ′ , each vertex has an edge (in H ′ ) with both its left neighbor and right neighbor in the arrangement.

Let

x 1 , x 2 ∈ V (H ′ ), where x 1 ∈ V (Q i ), x 2 ∈ V (Q j ).
We define the clockwise distance from x 1 to x 2 , denoted by clock dist(x 1 , x 2 ), as the minimum non-negative integer s such that j = i + s. Similarly, the clockwise distance from x 2 to x 1 , denoted by clock dist(x 2 , x 1 ), is the minimum non-negative integer s ′ such that i = j + s ′ . Let x, y, z ∈ V (H ′ ). We say y < x z if scanning the vertices of H ′ in clockwise direction along the circle starting from x, vertex y is encountered before z.

Let x ∈ V (Q i ). Vertex y is called the farthest neighbor of x before z if y ∈ N H ′ (x), y ∈ V (Q i ) ∪ V (Q i+1 ) ∪ V (Q i+2 ), y < x z, and for every other w ∈ N H ′ (x) either z < x w or w / ∈ V (Q i ) ∪ V (Q i+1 ) ∪ V (Q i+2 ) or both.
Claim 4.2. There always exists a vertex which is the farthest neighbor of x before z, unless

(x, z) ∈ E(H ′ ) and z ∈ V (Q i ) ∪ V (Q i+1 ) ∪ V (Q i+2 ).
(i) throughout this proof expressions involving subscripts of u, P , Q, and Q ′ are to be taken modulo p. Every such expression should be evaluated to a value in {0, . . . , p -1}. For example, consider a vertex u i , where i < p Then, p + i = i.

Let {A, B} be the bipartition of the bipartite graph H ′ . We categorize the edges of H ′ as follows: an edge Proof: Suppose H ′ has an l-edge, where l > 2, between x ∈ Q i and y ∈ Q i+l (see Fig. 1). Let

(x, y) ∈ E(H ′ ) is called an l-edge, if l = min(clock dist(x, y), clock dist(y, x)). b 2 a 1 b 1 a 2 u i u i-1 u i+l-1 u i+l y x Fig. 1: x ∈ V (Qi), y ∈ V (Q i+l )
a 1 = ||u i-1 Q ′ i x||, b 1 = ||xQ ′ i u i ||, a 2 = ||u i+l-1 Q ′ i+l y|| and b 2 = ||yQ ′ i+l u i+l ||.
We consider the following two cases: Case 1: l is even In this case u i-1 and u i+l-1 will be on the same side of the bipartite graph H ′ . Without loss of generality, let u i-1 , u i+l-1 ∈ A. Then, u i , u i+l ∈ B. We know that, for every

w 1 , w 2 ∈ V (H ′[m] ) with w 1 ∈ A and w 2 ∈ B, if (w 1 , w 2 ) / ∈ E(H ′[m] ) then d H ′ (w 1 , w 2 ) ≥ m + 2 (recalling m and d H ′ (w 1 , w 2 ) are odd integers). Therefore, we have a 1 + 1 + b 2 ≥ d H ′ (u i-1 , u i+l ) ≥ m + 2. Similarly, b 1 + 1 + a 2 ≥ d H ′ (u i , u i+l-1 ) ≥ m + 2. Summing up the two inequalities we get, (a 1 + b 1 ) + (a 2 + b 2 ) ≥ 2m + 2.
This implies that either ||Q ′ i || or ||Q ′ i+l || is greater than m which is a contradiction. Case 2: l is odd (proof is similar to the above case and hence omitted).

Hence we prove the claim.

We find a cycle C ′ = z 0 z 1 . . . z q z 0 in H ′ using Algorithm 3.1 (i) . Please read the algorithm before proceeding further. . (i) throughout this proof expressions involving subscripts of z are to be taken modulo q + 1. Every such expression should be evaluated to a value in {0, . . . , q}. For example, consider a vertex za, where a < q + 1. Then, q + 1 + a = a.

Algorithm 3.1 Finding Cycle C ′ in H ′ such that ||C ′ || ≥ ||C|| 1. l ← max l ′ (H ′ has an l ′ -edge).
Without loss of generality assume that this l-edge is between a vertex in Q 0 and a vertex in Q l 2. Scan the vertices of Q 0 in clockwise direction to find the first vertex z 0 , where z 0 ∈ V (Q 0 ), which has an l-edge to a vertex in Q l .

3. Scan the vertices of Q l in clockwise direction to find the last vertex in Q l which is a neighbor of z 0 in H ′ . Call it z 1 . 4. Find the farthest neighbor of z 1 before z 0 . Call it z 2 . /* refer proof of Claim 4.4 for a proof of existence of such a z 2 */ 5. s ← 2. while (z s , z 0 ) / ∈ E(H ′ ) do 6. Find the farthest neighbor of z s before z 0 . Call it z s+1 . /* such a neighbor exists by Claim 4.2*/ 7. s ← s + 1. end while 8. q ← s. 9. Return cycle C ′ = z 0 z 1 . . . z q z 0 . Claim 4.4. There always exists a farthest neighbor of z 1 before z 0 .

Proof: Note that z 0 ∈ Q 0 and z 1 ∈ Q l , where l ≤ 2 (by Claim 4.3). Recalling that ||C|| = p ≥ 6, we have z 0 / ∈ V (Q l ) ∪ V (Q l+1 ) ∪ V (Q l+2 ).
Hence by Claim 4.2, the claim is true.

Claim 4.5. The while loop in Algorithm 3.1 terminates after a finite number of iterations.

Proof: From Claim 4.1, we know that each vertex has an edge (in H ′ ) with both its left neighbor and right neighbor in the circular arrangement. Each time when Step 6 of Algorithm 3.1 is executed, a vertex z s+1 is chosen such that z s+1 is the farthest neighbor of z s before z 0 . Since H ′ is a finite graph, there will be a point of time in the execution of the algorithm when in Step 6 it picks a z s+1 such that (z s+1 , z 0 ) ∈ E(H ′ ).

From Claim 4.5, we can infer that C ′ is a cycle. Claim 4.6. C ′ is an induced cycle in H ′ .

Proof: Suppose C ′ is not an induced cycle. Then there exists a chord

(z a , z b ) in C ′ . Since (z a , z b ) is a chord, we have b = a -1 or b = a + 1. Let l = max l ′ (H ′ has an l ′ -edge). Let z a ∈ V (Q i ), z b ∈ V (Q j ). We know that min(clock dist(z a , z b ), clock dist(z b , z a )) ≤ l. Without loss of generality, assume clock dist(z a , z b ) ≤ l ≤ 2 (from Claim 4.3). That is, j -i ≤ l ≤ 2 and (z a , z b ) is a (j -i)-edge.
If z a = z 0 , then z b = z 1 and the algorithm exits from the while loop, when q = b, thus returning a cycle z 0 . . . z b z 0 . But in such a cycle (z b , z 0 ) is not a chord. Therefore, z a = z 0 . Similarly, z b = z 0 . We know that z a+1 = z b , z a+1 < za z b , and

z a+1 ∈ V (Q i ) ∪ V (Q i+1 ) ∪ V (Q i+2 ). Since j -i ≤ 2, z b ∈ V (Q i ) ∪ V (Q i+1 ) ∪ V (Q i+2 ). If z b < za z 0 ,
then it contradicts the fact that z a+1 is the farthest neighbor of z a before z 0 . Therefore, z 0 < za z b . Then, either z b = z 1 or z 1 < za z b . Recall that l = max l ′ (H ′ has an l ′ -edge), and (z 0 , z 1 ) is an l-edge with z 0 ∈ V (Q 0 ) and z 1 ∈ V (Q l ). Since (i) (z a , z b ) is a (ji)-edge, where ji ≤ l, (ii) z 0 < za z b , and (iii) (z a , z b ) is an l-edge. We know that (z 0 , z 1 ) is also an l-edge with

z b = z 1 or z 1 < za z b , we have l ≥ j -i = clock dist(z a , z b ) ≥ clock dist(z 0 , z b ) ≥ clock dist(z 0 , z 1 ) = l. Hence, j -i = l and
z 0 ∈ V (Q 0 ) and z 1 ∈ V (Q l ). Since z 0 < za z b and z b = z 1 or z 1 < z0 z b , we get z a ∈ V (Q 0 ) and z b ∈ V (Q l ). From
Step 2 of the algorithm we know that z 0 is the first vertex (in a clockwise scan) in Q 0 which has an l-edge to a vertex in Q l . This implies that, since z 0 < za z b , z a = z 0 which is a contradiction. Hence we prove the claim.

What is left now is to show that q + 1 ≥ p, i.e., ||C ′ || ≥ ||C||, where C ′ = z 0 . . . z q z 0 and C = u 0 . . . u p-1 u 0 . In order to show this, we state and prove the following claims. Claim 4.7. For every j ∈ {0, . . . , p -1}, (V (Q

j ) ∪ V (Q )) ∩ V (C ′ ) = ∅.
Proof: Suppose the claim is not true. Find the minimum j that violates the claim. Clearly, j = 0 as z 0 ∈ V (Q 0 ). We claim that z q ∈ V (Q j-1 ). Suppose z q / ∈ V (Q j-1 ). Let a = max{i | z i ∈ V (Q j-1 )} (note that, since j = 0, by the minimality of j, (V (Q j-1 )∪V (Q j ))∩V (C ′ ) = ∅ and therefore V (Q j-1 ) ∩ V (C ′ ) = ∅). Since z a = z q , by the maximality of a, we have

z a+1 / ∈ V (Q j-1 ). From our assumption, (V (Q j ) ∪ V (Q j+1 )) ∩ V (C ′ ) = ∅ and therefore z a+1 / ∈ V (Q j-1 ) ∪ V (Q j ) ∪ V (Q j+1
). Thus z a = z q and z a+1 is not the farthest neighbor of z a before z 0 . This is a contradiction to the way z a+1 is chosen by Algorithm 3.1. Hence, z q ∈ V (Q j-1 ). We know that (z q , z 0 ) ∈ E(H ′ ) with

z q ∈ V (Q j-1 ) and z 0 ∈ V (Q 0 ). Since l = max l ′ (H ′ has an l ′ -edge), we have min(clock dist(z q , z 0 ), clock dist(z 0 , z q )) ≤ l . That is, j ≥ p + 1 -l or j ≤ 1 + l. As l ≤ 2 (by Claim 4.3), we have j = p -1 or j ≤ 1 + l. Since z 0 ∈ V (Q 0 ), (V (Q p-1 ) ∪ V (Q 0 )) ∩ V (C ′ ) = ∅ and hence j = p -1. Therefore, j ≤ 1 + l. Since z 0 ∈ V (Q 0 ) and z 1 ∈ V (Q l ) (recall l ≤ 2), we get j = 1 + l. We know that, for every z a , z b ∈ V (C ′ ), if a < b then z a < z0 z b . Therefore, z 1 < z0 z q . We have z 1 ∈ V (Q l ). Since j = 1 + l, we also have z q ∈ V (Q l )
. Thus, we have z 1 , z q ∈ V (Q l ) and z 1 < z0 z q . But this contradicts the fact that z 1 is the last vertex in Q l encountered in a clockwise scan that has z 0 as its neighbor.

Claim 4.8. Let (z a , z a+1 ), (z b , z b+1 ) ∈ E(C ′ ) be two 2-edges, where a < b. Let P , P ′ denote the clockwise z a+1z b , z b+1z a paths respectively in C ′ . Both P and P ′ contain at least one 0-edge. Proof: Consider the path P (proof is similar in the case of path P ′ ). Path P is a non-trivial path only if z a+1 = z b . Suppose z a+1 = z b (see Fig. 2). Let z a ∈ V (Q f ). For the sake of ease of notation, assume f = 1 (the same proof works for any value of f ). Let

a 1 = ||u 0 Q ′ 1 z a ||, b 1 = ||z a Q ′ 1 u 1 ||, a 2 = ||u 2 Q ′ 3 z b ||, b 2 = ||z b Q ′ 3 u 3 ||, a 3 = ||u 4 Q ′ 5 z
c a+1 + c a+2 ≤ m -1 (1)
Since (u s+2 , u s+5 ) / ∈ E(H ′[m] ), (u s+1 , u s+4 ) / ∈ E(H ′[m]
), we have,

c b + d b+1 ≥ m + 1 c a+s + c b+1 ≥ m
Summing up the two inequalities, we get

c b + c a+s ≥ 2m + 1 -(c b+1 + d b+1 )
Since b = a + s + 1 and c b+1 + d b+1 ≤ m, we get

c a+s+1 + c a+s ≥ m + 1 (2)
Substituting for s = 1 in Inequality Claim 4.9. For every j, j ′ ∈ {0, . . . , p -1}, where j < j ′ and (

V (Q j ) ∪ V (Q j ′ )) ∩ V (C ′ ) = ∅, there exist i, i ′ ∈ {0, . . . , p -1}, where only i satisfies j < i < j ′ , such that |V (Q i ) ∩ V (C ′ )| ≥ 2 and |V (Q i ′ ) ∩ V (C ′ )| ≥ 2.
Proof: By Claim 4.7, (i) j ′ = j + 1 or j ′ = j -1, and (ii) there exist r, r ′ ∈ {0, . . . , q} such that (z r , z r+1 ) is a 2-edge with its endpoints on Q j-1 and Q j+1 and (z r ′ , z r ′ +1 ) is a 2-edge with its endpoints on Q j ′ -1 and Q j ′ +1 . By Claim 4.8, we know that if P , P ′ denote the clockwise z r+1z r ′ , z r ′ +1z r paths respectively in C ′ , then both P and P ′ contains at least one 0-edge. This proves the claim.

In order to show that the size of cycle C ′ (= z 0 . . . z q z 0 ) is at least p, we consider the following three cases:-Case |{Q j ∈ {Q 0 . . . Q p-1 } | V (Q j ) ∩ V (C ′ ) = ∅}| = 0: In this case, for every j ∈ {0, . . . p -1}, Q j contributes to V (C ′ ) and therefore ||C ′ || ≥ p = ||C||.

Case |{Q j ∈ {Q 0 . . . Q p-1 } | V (Q j ) ∩ V (C ′ ) = ∅}| = 1: Let Q j be that only path (among Q 0 . . . Q p-1 ) that does not contribute to V (C ′ ). Then we claim that there exists a Q j ′ , where j ′ = j, such that V (C ′ ) ∩ V (Q j ′ ) ≥ 2. Suppose the claim is not true then it is easy to see that ||C ′ || = p -1 which is an odd number thus contradicting the bipartitedness of H ′ . Hence the claim is true. Now, by applying the claim it is easy to see that

||C ′ || = j |V (C ′ ) ∩ V (Q j )| ≥ p = ||C||. Case |{Q j ∈ {Q 0 . . . Q p-1 } | V (Q j ) ∩ V (C ′ ) = ∅}| > 1:
Scan vertices of H ′ starting from any vertex in clockwise direction. Claim 4.9 ensures that between every Q j and Q j ′ , which do not contribute to V (C ′ ), encountered there exists a Q i which compensates by contributing at least two vertices to V (C ′ ). Therefore, ||C ′ || ≥ p = ||C||.

Discussion

An interesting open question that naturally follows from our result is the following: given a graph G and positive integers k, m where k ≥ 4, if G [m] is k-chordal, then is G [m+2] also k-chordal? As mentioned earlier, Brandstädt et al. in [START_REF] Brandstädt | Duchet-type theorems for powers of hhd-free graphs[END_REF] showed a similar result in the context of ordinary graph powering. They showed that, for every graph G, if G m is k-chordal, then so is G m+2 , where k, m are positive integers with k ≥ 3. A straightforward extension of their proof technique doesn't seem to work in our context due to the bipartite nature of the powering that we consider.

  Fig.1:x ∈ V (Qi), y ∈ V (Q i+l ) and let (x, y) ∈ E(H ′ ) be an l-edge, where l > 2. The dotted line between ui-1 and ui indicate the path Qi. Similarly, the dotted line between u i+l-1 and u i+l indicate the path Q i+l .

Fig. 2 :

 2 Fig. 2: Figure illustrates the case when path P defined in Claim 4.8 is a trivial path. The dotted lines between each ui-1 and ui indicate the path Q ′ i . Each continuous arc corresponds to an edge in the cycle C ′ = z0 . . . zqz0.

Fig. 3 :

 3 Fig. 3: Figure illustrates the case when path P defined in Claim 4.8 is P = za+1za+2 . . . za+1+s, where s ≥ 1 and za+1+s = z b . The dotted lines between each ui-1 and ui indicate the path Q ′ i . Each continuous arc corresponds to an edge in the cycle C ′ = z0 . . . zqz0.

  b+1 ||, and b 3 = ||z b+1 Q ′ 5 u 5 ||. We know that, for every w 1 , w 2 ∈ V (H ′[m] ) with w 1 ∈ A andw 2 ∈ B, if (w 1 , w 2 ) / ∈ E(H ′[m] ) then d H ′ (w 1 , w 2 ) ≥ m + 2. Since (u 0 , u 3 ) / ∈ E(H ′[m] ), (u 1 , u 4 ) / ∈ E(H ′[m] ) and (u 2 , u 5 ) / ∈ E(H ′[m] ), we have a 1 + b 2 ≥ m + 1, b 1 + a 3 ≥ m, and a 2 + b 3 ≥ m + 1.Adding the three inequalities and by applying an easy averaging argument we can infer that eithera 1 + b 1 = ||Q 1 || > m, a 2 + b 2 = ||Q 3 || > m, or a 3 + b 3 = ||Q 5 || > mwhich is a contradiction. Therefore P is a non-trivial path i.e., z a+1 = z b . Assume P does not contain any 0-edge. Let P = z a+1 z a+2 . . . z a+1+s , where s ≥ 1, a + 1 + s = b, and (z a+1 , z a+2 ) . . . (z a+s , z a+1+s ) are 1-edges (see Fig.3). Since (u0 , u 3 ) / ∈ E(H ′[m] ), (u 1 , u 4 ) / ∈ E(H ′[m] ),we have c a + d a+1 ≥ m + 1 and d a + d a+2 ≥ m (please refer Fig. 3 for knowing what c a , d a , . . . , c b+1 , d b+1 are). Summing up the two inequalities, we get d a+1 + d a+2 ≥ 2m + 1 -(c a + d a ). We know that, for each i ∈ {0, . . . p -1}, ||Q ′ i || ≤ m. Therefore, we have c a + d a ≤ m. Hence, d a+1 + d a+2 ≥ m + 1. Since (c a+1 + d a+1 ) + (c a+2 + d a+2 ) ≤ 2m, we get

  2, we get c a+2 + c a+1 ≥ m + 1. But this contradicts Inequality 1. Hence s > 1. Suppose s = 2. Since (u 2 , u 5 ) / ∈ E(H ′[m] )), we have c a+1 + d a+3 ≥ m. Adding this with Inequality 2, we get c a+1 + c a+2 ≥ (2m + 1) -(c a+3 + d a+3 ) ≥ m + 1. But this contradicts Inequality 1. Hence s > 2. Since (u s , u s+3 ) / ∈ E(H ′[m] )), . . . , (u 2 , u 5 ) / Adding the above set of inequalities and applying the fact that c i + d i ≤ m, ∀i ∈ {0, . . . q}, we get c a+1 + c a+2 + d a+s + d a+s+1 ≥ 2m. Adding this with Inequality 2, we get c a+1 + c

			∈ E(H ′[m] )), we have the
	following inequalities:-
	c a+s-1 + d a+s+1 ≥ m
	. . .	. . .	. . .
	c a+1 + d a+3 ≥ m

a+2 ≥ (3m + 1) -(c a+s+1 + d a+s+1 ) -(c a+s + d a+s ) ≥ m + 1. But this contradicts Inequality 1. Hence we prove the claim.