Claude Lobry

Tewfik Sari

Karim Yadi
email: k.yadi@mail.univ-tlemcen.dz

Coexistence of three predators competing for a single biotic resource

Keywords: 2000 Mathematics Subject Classification. 92A15, 92A17, 34C15, 34C35, 34E18 competition, coexistence, biotic resource, averaging, nonstandard analysis

We construct a model of competition of three consumers for one single biotic resource ; simulations show that the three species coexist. Using singular perturbations theory we sketch a mathematical proof for this coexistence. The main mathematical tool used is an extension of the Pontryagin-Rodygin theorem on the "slow" motion of a "slow-fast" differential system when the "fast" motion possesses a stable limit cycle. The mathematical analysis is done within the framework of Non Standard Analysis.

Introduction

The question of coexistence of competing species for a single resource has a very long history that we shall not attempt to recall here. We just recall the two decisive papers by Armstrong and Mac Gehee [START_REF] Armstrong | Coexistence of species competing for shared resources[END_REF][START_REF] Mcgehee | Some mathematical problems concerning the ecological principle of competitive exclusion[END_REF] where they pointed that coexistence is not synonymous of coexistence at equilibrium. These papers were the starting point of numerous papers showing complex behaviors of systems of competitors and evidence of coexistence on the basis of numerical simulations. Following this tradition we propose a model of coexistence of three species competing for one resource.

The present paper has two parts. In the first part we construct our model and explain what is the rationale behind our construction ; each step is illustrated by simulations. In the second part we consider our model as a member of a more general "consumer-resource" model for which we explain how coexistence of species can be proved using singular perturbation analysis ; as an essential tool we use an extension of a theorem of Pontryagin and Rodygin.

Considered from the ecological point of view our model shows that oscillations in a "consumer-resource" relationship can open the door to coexistence with other species provided that the new introduced species do not perturb too much the oscillations. We do not know any example of an interaction between four species of the type of the model presented here but its existence is plausible. We shall explain it during the construction of a model. But we must acknowledge here that what we do is a kind of "virtual ecology" showing what is "theoretically possible" in a world of species respecting basic facts well established in concrete ecology. It is not a description of the real world !

From the mathematical point of view our paper can be considered as an application of singular perturbation methods to the mathematical proof of persistence for some specific system. Our contribution consists mainly in the analysis of the theorem of Pontryagin and Rodygin and its extension to a theorem which is more effective in some circumstances. Detailed proofs can be consulted at [START_REF] Yadi | Perturbations Singulières : Approximations, Stabilité Pratique et Applications à des Modèles de Compétition[END_REF] and will be published elsewhere. The mathematical analysis is done within the framework of Non Standard Analysis, using the axiomatic of Nelson [START_REF] Nelson | Internal Set Theory: a new approach to nonstandard analysis[END_REF] and respecting the spirit of G. Reeb [START_REF] Reeb | La mathématique non standard vieille de soixante ans ? Troisième Colloque sur les Catégories, dédié à C. Ehresmann, Amiens[END_REF].

2. Construction of a model.

2.1.

The basic oscillating pair. We consider the system :

(1)

     ds dt = 3s(1 -s) -s 2 0.01+s 2 x 1 dx 1 dt = 0.1[ s 2 0.01+s 2 -0.65]x 1
Due to the presence of the factor 0.1, it can be considered as a "slow-fast" system of two differential equations of the following type :

     ds dt = 1 ε1 [f (s) -g 1 (s)x 1 ] dx 1 dt = (g 1 (s) -d 1 )x 1
where ε 1 is a "small" parameter. The real s represents the density of some biotic resource (prey) for a consumer which density is represented by x 1 . This is a rather classical prey-predator model and it is well known that this kind of generalization of the Lotka-Volterra model can have sustained oscillations. The existence of oscillations in prey-predator interaction is clearly demonstrated in laboratory experiments (see [START_REF] Jost | Interaction of Tetrahimen pyriformis, Escherichi coli, Azotobacter vinelandi and Glucose in a Minimal Madium[END_REF]), if not in the real world. The choice of the function :

g 1 (s) = s 2 0.01 + s 2
in place of the more classical Monod's function with s in place of s 2 was made in order that the nullcline [f (s) -g 1 (s)x 1 ] = 0 has the "S-shape" shown on Fig. 1 which prevents the resource from extinction. This kind of assumption is sometimes called the "Allee" effect in ecological literature.

2.2. Addition of a new consumer "x 2 ". We want to add a new consumer and at the same time keep the oscillations of (s, x 1 ). It can be done by introducing a new species with a very slow dynamics compare to that of (s, x 1 ) like in the model :

             ds dt = 1 ε1 [f (s) -g 1 (s)x 1 -g 2 (s)x 2 ] dx 1 dt = (g 1 (s) -d 1 )x 1 dx 2 dt = ε 2 (g 2 (s) -d 2 )x 2
where ε 2 is small. The existence of two consumers (which densities are represented by x 1 and x 2 ) of the same resource and having very different characteristic time seems to be common in nature. For instance small mammals and big mammals eating the same grass have a lifespan which may differ of an order of magnitude. Denote by x 2 some constant x 2 ; since x 2 (t) is quasi constant, for a while, the evolution of (s, x 1 ) is governed by :

       ds dt = 1 ε1 [f (s) -g 1 (s)x 1 -g 2 (s)x 2 ] dx 1 dt = (g 1 (s) -d 1 )x 1
In that system we see that when x 2 is small the nullcline

(2) f (s) -g 1 (s)x 1 -g 2 (s)x 2 = 0 is very close to the nullcline (3) f (s) -g 1 (s)x 1 = 0
and, thus, oscillations are preserved. The range of oscillations of s is slightly shortened as x 2 increases. Now, if we look at the process from the point of view of x 2 during an oscillation of period T the growth is given by :

t+T t (g 2 (s(τ )) -d 2 )x 2 (τ )dτ
which varies monotonically according to the variation of amplitude of s. Thus the growth can be positive for small values of x 2 and negative for large ones ; in the Figure 2. Simulation of the system (4). The picture shows the two projections of the trajectories, in green on the plane (s, x 1 ) and in red on the plane (s, x 2 ). While the initial conditions where kept constant for s and x 1 they where changed for x 2 . One sees that, starting from a small x 2 = a then x 2 (t) is increasing and, conversely, starting from a big x 2 = c then x 2 (t) is decreasing.

middle there must be an equilibrium. This is the case for the model :

(4)              ds dt = 3s(1 -s) -s 2 0.01+s 2 x 1 -g 2 (s)x 2 dx 1 dt = 0.1[ s 2 0.01+s 2 -0.65]x 1 dx 2 dt = 0.01[g 2 (s) -0.025]x 2 with s < 0.5 ⇒ g 2 (s) = 0 s ≥ 0.5 ⇒ g 2 (s) = 0.1(s-0.5) 0.01+(s-0.5)
We choose g 2 = 0 on [0, 0.5] in order to be sure that the nullcline (2) is exactly the same than the nullcline (3) and remains "S-shaped" ; thus oscillations are preserved. This artificial choice is convenient for simplicity but a model with g 2 being smoother would also work. Evidence of coexistence of x 1 and x 2 is given on Fig. 2. On this simulation three initial conditions where taken, keeping s(0) and x 1 (0) constant and changing x 2 (0). On the picture we have superimposed the projections on the (s, x 1 ) plane (in green) and the (s, x 2 ) plane in red. The first initial condition for x 3 is a which is shown by the black arrow ; since the variation of x 2 is very slow the red projection looks like a point moving right to left and left to right very fast while moving up very slowly ; the thick red line corresponds to a dozen of oscillations. Starting from c we see that the point is slowly moving down.

Starting from b we see that x 2 remains constant.

2.3. Addition of a new consumer "x 3 ". The idea is to have a new consumer with a g 3 growth rate such that the two graphs of g 2 and g 3 cross like on Fig.

3 in order that during an oscillation of (s, x 1 ) the species 2 and 3 take the advantage alternatively. This leads to the system :

(5)

             ds dt = 3s(1 -s) -s 2 0.01+s 2 x 1 -g 3 (s)x 3 dx 1 dt = 0.1[ s 2 0.01+s 2 -0.65]x 1 dx 3 dt = 0.01[g 3 (s) -0.025]x 3 with s < 0.58 ⇒ g 3 (s) = 0 s ≥ 0.58 ⇒ g 3 (s) = 2(s-0.58) 0.01+(s-0.58)
On Fig. 4 one sees that the behavior of the system (s, x 1 , x 3 ) is similar to the behavior we observed for (s, x 1 , x 2 ) ; the projection (in blue) on the (s, x 3 ) plane is similar to the projection (in red) observed for (s, x 2 ) in Fig. 2.

2.4.

Coexistence of all the three consumers. Now we consider the complete system :

             ds dt = 3s(1 -s) -s 2 0.01+s 2 x 1 -g 2 (s)x 2 -g 3 (s)x 3 dx1 dt = 0.1[ s 2 0.01+s 2 -0.65]x 1 dx2 dt = 0.01[g 2 (s) -0.025]x 2 dx3 dt = 0.01[g 3 (s) -0.025]x 3 (6)
Species x 2 and x 3 have a "slow motion" which can be approximated by computing suitable integrals along the basic cycle ; this determines a flow on (x 2 , x 3 ) plane ; this flow is studied and proved to have a stable equilibrium which proves the persistence of both x 2 and x 3 . On Fig. 5 one sees the projection on the (x 2 , x 3 ) plane of the full system (6) with the three competitors. The simulations from various initial conditions shows convergence to a point which actually corresponds to a periodic orbit in the full space (s, x 1 , x 2 , x 3 ). More details are given in the next section.

2.5. Species 2 and 3 alone. Consider the system with s, x 2 and x 3 alone in the absence of the species represented by x 1 . One easily check that in this case, from the choice of g 2 and g 3 , there is no oscillation and the stable equilibrium is the one for which species x 2 wins the competition. 5). The picture shows the two projections of the trajectories, in green on the plane (s, x 1 ) and in blue on the plane (s, x 3 ). The projection (in blue) on the (s, x 3 ) plane is similar to the projection (in red) observed for (s, x 2 ) in Fig. 2. Starting from a small x 3 = a then x 3 (t) is increasing and, conversely, starting from a big x 3 = c then x 3 (t) is decreasing. Between a and c there is some initial condition (not represented) for which x 3 do not increase nor decrease.

Figure 5. On this picture we have represented the projection of simulations of the the complete system (6) on the plane (x 2 , x 3 ). The variables (s 1 , x 1 ) (not represented) present rapid oscillations while (x 2 , x 3 ) evolves slowly. One sees that all the trajectories seem to converge to an equilibrium. Compare to the "theoretical"picture on Fig 8.

Proof of the persistence of the three competitors in a model with three time scales

In this section, we give the successive steps of a proof of the persistence in a model of the form ( 7)

       ds dt = 1 ε1ε2 (f (s) -g 1 (s)x 1 -g 2 (s)x 2 -g 3 (s)x 3 ) dx1 dt = 1 ε1 (g 1 (s) -d 1 )x 1 dx2 dt = (g 2 (s) -d 2 )x 2 dx3 dt = (g 3 (s) -d 3 )x 3
where the occurring functions are differentiable, at least piecewise, ε 1 and ε 2 are infinitesimal 1 . The function f vanishes at 0 ; it is increasing then decreasing and vanishes at a value m. The functions g i are zero at 0, increasing and bounded.

3.1. Oscillations of s and x 1 . Let us consider the system ( 8)

ds dt = 1 ε1 (f (s) -g 1 (s)x 1 ) dx1 dt = (g 1 (s) -d 1 )x 1
• Suppose that the nullcline ds/dt = 0 is a curve ϕ that, when s increases from 0, decreases from +∞ to a minimum value reached for s = s -then increases to a maximum value for s = s + to finally decrease and vanishes for s = m. • Suppose that the value s * such that g 1 (s * ) = d 1 is between s -and s + . Proposition 3.1. For ε 1 infinitesimal, the system (8) has a limit cycle close to the curve ABCD in Fig. 6.

3.2.

The Pontryagin-Rodygin's Theorem. Due to the lack of space, we shall try in these sections to avoid excessive mathematical formalism of the results and we refer to [START_REF] Sari | On Pontryagin-Rodygin's theorem for convergence of solutions of slow and fast systems[END_REF][START_REF] Yadi | Perturbations Singulières : Approximations, Stabilité Pratique et Applications à des Modèles de Compétition[END_REF][START_REF] Yadi | Singular perturbations on infinite time interval[END_REF] for more details. To simplify, we suppose that all the occurring differential equations have the property of uniqueness of solutions. Let us consider the slow and fast system

(9)      ds dt dx1 dt = 1 ε F (s, x 1 , x) G(s, x 1 , x) dx dt = H(s, x 1 , x)
where the scalars s and x 1 are the fast components, and the vector x the slow one.

The real number ε is positive and infinitesimal. The functions F , G and H are continuous. The following system, where x is considered as a parameter, is called the fast equation

ds dt dx1 dt = 1 ε F (s, x 1 , x) G(s, x 1 , x) (10)
Hence, the (s, x 1 )-component of a solution of (9) varies very quickly according to [START_REF] Robinson | Nonstandard Analysis[END_REF] where x has been frozen at its initial value. When the fast equation ( 10) has stable limit cycles Γ x for each x in a compact domain, Pontryagin-Rodygin's 1 We use N.S.A. terminology. Following the spirit of ( [START_REF] Lobry | Nonstandard analysis and representation of real world[END_REF]) when a word (like "infinitesimal") is written in bold character its meaning is the one used in the formal language of Nelson I.S.T. but the reader not familiar with this framework can use the intuitive meaning.

G=0 F=0 s=! 1 (x 1 ,x) s=! 2 (x 1 ,x) "1 (x) " 2 (x) s' >0 x'1>0 s' >0 x'1<0 s' <0 x'1>0 s' <0 x'1<0 s x 1 Figure 7.
Notations in equations ( 14) admits a stable limit cycle for any value of x for infinitesimal values of ε 1 . It is tempting to apply Pontryagin-Rodygin's Theorem to [START_REF] Sari | On Pontryagin-Rodygin's theorem for convergence of solutions of slow and fast systems[END_REF] but the main reason that makes this impossible is the fact that the fast equation is a nonstandard equation. It is itself a singularly perturbed equation. We can not avoid to take into account the three dynamics of the problem. In [START_REF] Yadi | Perturbations Singulières : Approximations, Stabilité Pratique et Applications à des Modèles de Compétition[END_REF] Pontryagin-Rodygin's Theorem is extended to this kind of system the fast equation of which admits a slow and fast limit cycle. This new result has also the advantage to overcome a serious limitation of Pontryagin-Rodygin's Theorem : unlike the latter, it makes possible the localization of the cycles, the approximation of their periods and the calculation of the average along these cycles. The functions F, G and H being continuous and the positive real numbers ε 1 and ε 2 infinitesimal, suppose that there exists a compact domain K of R 2 such that, for all x ∈ K the nullclines F = 0 and G = 0 of (13) have the shape given in Figure 7. The (s, x 1 )-plane is divided in four regions where the field has the indicated signs in the figure. The limit cycle of ( 13) is infinitesimally close to the closed curve (ABCD) in Fig. 6 formed by two "slow arcs" (AB) and (CD) and two "fast segments" (DA) and (BC). The two decreasing branches of the nullcline F = 0 are denoted s = ψ 1 (x 1 , x) and s = ψ 2 (x 1 , x). Let us define in the interior of K the slow equation

(14) dx dt = M (x), where (15)
M (x) = 1 P (x) 2 i=1 ξi+1(x) ξi(x) g(ψi(x1,x),x1,x) f2(ψi(x1,x),x1,x) dx 1 , P (x) = 2 i=1 ξi+1(x) ξi(x) dx1 f2(ψi(x1,x),x1,x) , with ξ 3 (x) = ξ 1 (x).
Let γ(t) be the trajectory of a solution of (12). Theorem 5.2.1 page 75 in [START_REF] Yadi | Perturbations Singulières : Approximations, Stabilité Pratique et Applications à des Modèles de Compétition[END_REF] • Assume that the functions g 2 and g 3 are zero until the respective thresholds s 2 and s 3 are reached such that min(s 2 , s 3 ) ≥ s + and that g 2 and g 3 are increasing beyond. This assumption allows us to assert that the subsystem

ds dt dx1 dt = 1 ε 2 1 ε1 (f (s) -g 1 (s)x 1 -g 2 (s)x 2 -g 3 (s)x 3 ) (g 1 (s) -d 1 )
x 1 still admits, for every (x 2 , x 3 ) and ε 1 small enough, a limit cycle Γ x2,x3 that differs from that of (8) for values of s ≥ min(s 2 , s 3 ). The more x 1 and x 2 are large, the more these cycles are distorted inwards in their right side. The minimum and maximum of the cycles remain unchanged. Here, the averaged equation [START_REF] Yadi | Perturbations Singulières : Approximations, Stabilité Pratique et Applications à des Modèles de Compétition[END_REF]15) takes the form (see [START_REF] Yadi | Perturbations Singulières : Approximations, Stabilité Pratique et Applications à des Modèles de Compétition[END_REF] for explicit formulas) (16) dx 2 /dt = x 2 M 2 (x 2 , x 3 )/P (x 2 , x 3 ), dx 3 /dt = x 3 M 3 (x 2 , x 3 )/P (x 2 , x 3 ).

A detailed study of equation ( 16) leads to the following conditions of persistence : Theorem 3.2. [START_REF] Yadi | Perturbations Singulières : Approximations, Stabilité Pratique et Applications à des Modèles de Compétition[END_REF] Suppose that s 3 > s 2 = s + , ϕ(s 3 ) > ϕ(s -) and that d 2 is below a certain constant well determined by the problem. Then, for s 3 -s 2 and d 3 small enough, there is persistence of the species x 2 and x 3 of (16).

This result is reflected in Fig. 8 representing a positively invariant box OABC of (16) in which arrive all trajectories with positive initial conditions. The axes are invariant, the origin O is an unstable node and C and E are saddle points. A lemma due to Butler-McGehee [START_REF] Butler | Persistence in dynamical systems[END_REF] shows that the union of limit sets of positive half-trajectories is a compact subset Ω which does not meet the axes. Theorem 3.3. Under the assumptions of the preceding Theorem, there is persistence of the whole species of the model [START_REF] Nelson | Internal Set Theory: a new approach to nonstandard analysis[END_REF] for all positive initial conditions. This final result is obtained by using a nonstandard permanence lemma [START_REF] Robinson | Nonstandard Analysis[END_REF] which extends the approximation of the component (x 2 (t), x 3 (t)) of the solution of (7) by the solution of (16) to an infinitely large time interval [0, ω]. We than prove2 that (x 2 (t), x 3 (t)) remains infinitely close to Ω for all infinitely large values of t and all infinitesimal values of ε 1 and ε 2 .

Figure 1 .

 1 Figure1. On this picture one observes a simulation of few trajectories of the system (1). The black "S-shaped"curve is the nullcline [f (s) -g 1 (s)x 1 ] = 0. The direction of the motion along trajectories is indicated by the black arrows and the limit cycle by the red arrow. By the way s(t) oscillates between two values.

Figure 3 .

 3 Figure 3. The graphs of g 2 and g 3 .

Figure 4 .

 4 Figure 4. Simulation of the system[START_REF] Lobry | Nonstandard analysis and representation of real world[END_REF]. The picture shows the two projections of the trajectories, in green on the plane (s, x 1 ) and in blue on the plane (s, x 3 ). The projection (in blue) on the (s, x 3 ) plane is similar to the projection (in red) observed for (s, x 2 ) in Fig.2. Starting from a small x 3 = a then x 3 (t) is increasing and, conversely, starting from a big x 3 = c then x 3 (t) is decreasing. Between a and c there is some initial condition (not represented) for which x 3 do not increase nor decrease.

Figure 8 .

 8 Figure 8. Portrait of (x 2 , x 3 ) obtained from the averaged system (16). Compare to the simulations presented on Fig.5

3. 4 .

 4 Application to the model. Reconsider System[START_REF] Nelson | Internal Set Theory: a new approach to nonstandard analysis[END_REF].

We say that the standard set {Γx 2 ,x

× {(x 2 , x 3 )} : (x 2 , x 3 ) ∈ Ω} is practically asymptotically stable for[START_REF] Nelson | Internal Set Theory: a new approach to nonstandard analysis[END_REF] for all infinitesimal values of ε 1 and ε 2 (see[START_REF] Yadi | Singular perturbations on infinite time interval[END_REF] for more details).

Starting from the point P the solution is quasi horizontal and goes fast to the nullcline ds/dt = 0 ; then the solution stays near the nullcline and goes up until it reaches the maximum at B; then the solution is quasi horizontal and goes fast to the nullcline ds/dt = 0 at point C ; then the solution stays near the nullcline and goes down until it reaches the minimum at D ; then the solution goes fast to the right to the nullcline ds/dt = 0 and reaches it at A.

Theorem [START_REF] Pontryagin | Approximate solution of a system of ordinary differential equations involving a small parameter in the derivatives[END_REF] gives the limiting behavior of the singularly perturbed problem (9) : Under suitable conditions, after a fast transition near the cycles described by the fast equation [START_REF] Robinson | Nonstandard Analysis[END_REF], the trajectories of ( 9) quickly roll up around the manifold generated by the cycles, with a slow evolution of the x-component according to the averaged system [START_REF] Sari | Averaging in Ordinary Differential Equations and Functional Differential Equations, The Strength of Nonstandard Analysis[END_REF] dx dt = 1 P (x)

where (s * (τ, x), x * 1 (τ, x)) is a P (x)-periodic solution of the fast equation corresponding to the cycle Γ x . This result was originally obtained for at least C 2 vector fields, under the assumption that the cycles Γ x are asymptotically stable in the linear approximation. However, the result obtained in [START_REF] Sari | On Pontryagin-Rodygin's theorem for convergence of solutions of slow and fast systems[END_REF] shows that Pontryagin-Rodygin description of solutions holds for C 0 vector fields under additional assumptions.

Extension of


Dr Julien Chevallier
email: j.chevallier@imperial.ac.uk

Modelling the convenience yield in carbon prices using daily and realized measures

Keywords: JEL Classification: C5, G1, Q4 Convenience Yield, Carbon Price, EU ETS, High-frequency Data, Realized Volatility

This article investigates the modelling of the convenience yield in the European carbon market by using daily and intradaily measures of volatility. The convenience yield stems from differences in spot and futures prices, and can explain why firms hold inventories. The main findings are that (i) a simple AR(4) process best describes the 2008 convenience yield, and (ii) there exists a non-linear relation between spot and futures prices. The approach developed in this article captures 74% of the explanatory power for the 2008 convenience yield variable in an autoregressive framework, with carbon spot price levels, moving averages and carbon futures realized volatility measures as exogenous regressors. These results are of interest for energy utilities, risk-managers, and traders exposed to the variation of carbon prices.

Introduction

On commodity markets, the investigation of the convenience yield appears as a central empirical issue, since it allows practioners, hedgers, brokers and other market operators in the field to buffer themselves against unanticipated changes in market conditions (see [START_REF] Pindyck | The dynamics of commodity spot and futures markets: a primer[END_REF] for a through discussion on this issue). The basic intuition is that differences in spot and futures prices occur due to the cost of holding inventory.

In many economic and financial applications, the theory of storage aims at explaining the differences between spot and futures prices by analyzing the reasons why agents hold inventories. According to [START_REF] Geman | Commodity and Commodity Derivatives: Modeling and Pricing for Agriculturals, Metals and Energy[END_REF], inventories have a productive value since they allow to meet unexpected demand, avoid the cost of frequent revisions in the production schedule and eliminate manufacturing disruption. [START_REF] Working | The theory of the price of storage[END_REF] defined the notion of convenience yield as a benefit that accrues to the owner of the physical commodity. [START_REF] Brennan | The supply of storage[END_REF] further defined the convenience yield as an embedded timing option attached to the commodity, since inventory allows to put the commodity on the market when prices are high and hold it when prices are low. Recent applications of the theory of storage to the modelling of commodity prices include [START_REF] Considine | Risk premium on inventory assets: the case of crude oil and natural gas[END_REF], [START_REF] Wei | Commodity convenience yield and risk premium determination: the case of the U.S. natural gas market[END_REF], [START_REF] Geman | Forward curves, scarcity and price volatility in oil and natural gas markets[END_REF] and [START_REF] Stronzik | Does the European natural gas market pass the competitive benchmark of the theory of storage? Indirect tests for three major trading points[END_REF] for crude oil or natural gas markets. This article focuses on the modelling of the convenience yield for carbon spot and futures prices, which are exchanged since 2005 on the European Union Emissions Trading Scheme (EU ETS). The EU emissions trading system has been created by the Directive 2003/87/CE. Across 27 Member States, the EU ETS covers large plants from CO 2 -intensive emitting industrial sectors with a rated thermal input exceeding 20 MWh.

One allowance exchanged on the EU ETS corresponds to one ton of CO 2 released in the atmosphere, and is called a European Union Allowance (EUA) (see [START_REF] Alberola | Price drivers and structural breaks in European carbon prices 2005-2007[END_REF] for more details.). This issue is of particular importance for risk-managers and traders in energy utilities regulated by the scheme, as they need to cover themselves against financial, political, and economic risks specific to this market (see Chevallier et al. (2009) for more details). Besides, carbon prices convey some interesting characteristics in terms of commodity modelling, since the costs of storage are null. Carbon permits indeed only exist in the balance sheets of the companies regulated by the scheme, and the costs of storing them physically are insignifiant.

There is very limited literature on the investigation of the convenience yield in the European carbon market. To our best knowledge, only [START_REF] Borak | Convenience Yields for CO 2 Emission Allowance Futures Contracts[END_REF] address this issue. They show that the market has changed from initial backwardation to contango with significant convenience yields in future contracts for the Kyoto commitment period starting in 2008. Their main result features that a high fraction of the yields can be explained by the price level and volatility of the spot prices. The authors conclude that the yields can be interpreted as market expectation on the price risk of CO2 emissions allowance prices and the uncertainty of EU allocation plans for the Kyoto period.

Based on recent developments in financial econometrics (Andersen, Bollerslev, Diebold and Labys, henceforth ABDL (2003)), we use in this article both daily and intradaily data for risk measures to model the convenience yield. High-frequency data is indeed superior in estimation if some biases (linked to microstructure noise) are correctly accounted for2 . This methodology, joint with standard ARMA filtering, yields some interesting results for the modelling of the convenience yield in carbon spot and futures prices, which is predictable (with a R-squared of 74%) using autoregressive processes. Using realized volatility measures shows two effects: 1) the explanatory power is reduced to 35%, and 2) realized volatility significantly impacts the convenience yield in carbon prices. Our results extend [START_REF] Borak | Convenience Yields for CO 2 Emission Allowance Futures Contracts[END_REF], who used daily data only and whose study period covers the early years of the EU ETS (2005)(2006).

The remainder of the article is composed as follows. Section 2 describes the data used. Section 3 details the modelling of the convenience yield. Section 4 develops realized volatility estimation techniques. Section 5 presents the empirical results. Section 6 concludes.

Data

We detail below the data used for carbon prices at both daily and intra-daily frequencies, as well as for the risk-free rate. [START_REF] Alberola | European carbon prices and banking restrictions: Evidence from Phase I (2005-2007)[END_REF], banking restrictions between 2007 and 2008 caused the disconnection between spot and futures prices between the two Phases. Besides, a structural break due to information revelation occurred in April 2006 for carbon prices of all maturities [START_REF] Alberola | Price drivers and structural breaks in European carbon prices 2005-2007[END_REF].).

Second, concerning intraday carbon futures prices, our sample contains one year of tick-by-tick transactions for the ECX futures contract of maturity December 2008, going from January 2 to December 15, 2008. This is equivalent to 240 days of trading after cleaning the data for outliers, and until the expiration of the contract. The average amount of transactions for the ECX carbon futures tick-data is equal to 700 trades per day. This corresponds to an average of 50 seconds between each transaction.

Third, the risk-free rate used below to compute the convenience yield between carbon spot and futures prices is the Euribor, as commonly used by market agents. The Euribor rates were accessed from Thomson Financial DataStream. Depending on the time until maturity, we use the Euribor contract with the relevant maturity. Descriptive statistics for all variables used in our econometric specification are given in Table 1.

We observe that CO 2 spot and futures price series are characterized by a negative skewness, and the kurtosis coefficient is close to three, which is the value for the normal distribution. Taken together, these descriptive statistics suggest that CO 2 spot and futures price series exhibit some leptokurticity, which can be better fitted by the used of GARCH(p, q) modeling [START_REF] Bollerslev | Generalized autoregressive conditional heteroskedasticity[END_REF]).

Next, we provide some elements on the theory of storage and develop the computational steps in order to obtain the convenience yield from carbon spot and futures prices.

Modelling the convenience yield

This section describes how the convenience yield can be measured.

By the cost-of-carry relationship, and without storage costs for EUA allowances, the futures and spot prices are linked through S t = F T e -r(T -t) with S t the spot price at time t, F T the futures prices of a contract with delivery in T and r the interest rate [START_REF] Working | The theory of the price of storage[END_REF], [START_REF] Brennan | The supply of storage[END_REF]). Equivalently, we may write:

F t,T = S t e (r-y)(T -t) (1)
with F t,T the futures price for maturity T at time t, r the continuously compounded risk-free interest rate used by market agents at time t for maturity T , and y is the convenience yield on the commodity.

As [START_REF] Pindyck | The dynamics of commodity spot and futures markets: a primer[END_REF] put it, this no-arbitrage condition states that the only cost of buying a commodity at time t and delivering it at maturity T is the foregone interest. Agents incur the opportunity cost of purchasing the asset, but in return they benefit from possessing the commodity and being able to trade it until maturity.

Hence, the convenience yield at time t for maturity T may be modelled directly as:

y t,T = S t e r(T -t) -F t,T (2)
Figure 2 shows the carbon spot and futures prices during 2008, as well as the corresponding convenience yield3 .

In the bottom panel of Figure 2, we may observe that the convenience yield is strongly time-varying during 2008, going from -0.2 in June to 1.2 in April. These variations may be explained by (i) the delayed effect of the "credit crunch" crisis on the carbon market [START_REF] Chevallier | Carbon futures and macroeconomic risk factors: A view from the EU ETS[END_REF]), and (ii) the 2007 compliance event which occurred in April 2008 (see Chevallier et al. (2009) for a detailed analysis of the effects of compliance events on investors' expectation changes).

To test this relationship empirically, we adopt the following econometric specification4 :

y t,T = α + β S t + γV t + ε t (3)
with V t a proxy for volatility, ε t the error term. y t,T is filtered through an ARMA process, which does not appear in the specification below, following the Box-Jenkins methodology. V t may be either an autoregressive or a moving average process for a measure of volatility, which may be composed of realized volatility measures (as in ABDL ( 2003)). Thus, we will explicitly compare in Tables 7 and8 several volatility proxies, which may be based on daily data (through the use of moving averages) or based on intraday data (through the use of realized volatility estimates). Since realized volatility estimates are computed over the time interval of one day, our econometric model may be carried out with a daily frequency.

In the next section, we detail how to compute realized volatility measures for carbon prices.

Estimation of realized volatility

This section presents different computation methods for realized volatility measures using the most recent developments in financial econometrics techniques.

Let p(t) denote a logarithmic asset price at time t. With no jump, the continuous-time diffusion process generally employed in asset and derivatives pricing may be expressed by a stochastic differential equation as:

d p(t) = µ(t)dt + σ (t)dW (t) with 0 ≤ t ≤ T (4)
with µ(t) a continuous and locally bounded variation process, σ (t) a strictly positive càdlàg (right continuous with left limits) stochastic volatility process, and W (t) a standard Brownian motion.

Next, let us consider the quadratic variation (QV) for the cumulative return process r(t) ≡ p(t)p(0):

[r, r] t = t 0 σ 2 (s)ds (5)
The QV simply equals the integrated volatility of the process described in Eq. ( 4). The realized volatility (RV) is defined as the sum of returns at a frequency 1/∆, or:

RV t+1 (∆) ≡ 1/∆ ∑ j=1 r 2 t+ j.∆,∆ (6)
When ∆ → 0, using theory of quadratic variations, it can be shown (see [START_REF] Andersen | The distribution of stock return volatility[END_REF], [START_REF] Barndorff-Nielsen | Econometric analysis of realized volatility and its use in estimating stochastic volatility models[END_REF]) that:

RV t+1 (∆) → t 0 σ 2 (s)ds (7)
Theory suggests that optimal sampling corresponds to sampling at the highest possible frequency. However, this is not true in practice due to microstructure effects (bid-ask spread, rounding, non-synchronicity, etc.) which introduce noise in the price process. To mitigate the impact of microstructure noise, we examine the volatility signature plot (as in ABDL ( 2003)) for the 2008 carbon futures contract. In volatility signature plots, the realized volatility measure described in Eq. ( 6) is computed and plotted at different sampling frequencies.

Figure 3 shows that as frequency becomes higher, the realized variance includes an increasing noise component. Thus, the optimal sampling frequency for carbon prices may be determined at 15-min returns (see [START_REF] Chevallier | Carbon futures and macroeconomic risk factors: A view from the EU ETS[END_REF]Sévi (2009a, 2009b) for a detailed analysis).

Figure 4 plots the three proxies of volatilities using realized measures. We use the last tick method which is superior to interpolation (ABDL ( 2003)). Subsampling [START_REF] Zhang | A tale of two time scales: determining integrated volatility with noisy high frequency data[END_REF]) and kernel-based (Zhou (1996), [START_REF] Hansen | Realized variance and market microstructure noise[END_REF]) estimators are similar in nature. The time-series reveal the presence of jumps and structural breaks that may be taken into account using bipower and tripower variation measures (ABD ( 2007)). The extension of our work to jump-robust measures of realized volatility is also left for further research.

In the next section, we present the estimation results of Eq. ( 3) for the 2008 convenience yield between carbon spot and futures prices using daily and realized measures.

Empirical results

Our empirical approach is of interest for policy makers, since it will allow to derive informed conclusions on the modelling of the convenience yield in the carbon market, which can then be used for forecasting purposes. Thus, we estimate eq(3) with various specifications.

Using daily and intraday data as a risk measure, we aim at identifying here whether the convenience yield in carbon prices is highly time-variant. Within this framework, volatility can either be an autoregressive or moving average process which can possibly consist of realized volatility measures.

In this section, we present the modelling results of Eq. ( 3) following different specifications of the 2008 convenience yield variable through ARMA filtering, and the inclusion as exogenous regressors of carbon spot price levels, moving averages, and carbon futures realized volatility estimates.

We first run the Augmented Dickey-Fuller (1979) unit root test on the dependent variable.

From Table 2, we may reject the unit root hypothesis for the 2008 convenience yield variable. The time-series does not seem to be integrated of any order. gressive components up to order 4 (AR( 4)): all lag orders are statistically significant at the 1% level. The R 2 is equal to 74%, which shows the high degree of predictability of the 2008 convenience yield variable for carbon spot and futures prices based on autoregressive processes only. All diagnostic tests are validated for this regression, and residuals are not autocorrelated5 .

Variance estimation is a critical part of this article. Thus, we need to account for heteroskedasticity, as if often the case for financial time-series. To detect ARCH-effects, we re-estimate this regression with an ARCH term in the variance equation, as shown in Table 4.

The ARCH component is significant at the 5% level, while other estimated coefficients remain stable.

This regression confirms the robustness of our previous results with an AR(4) process. This specification also allows to adjust for serial correlation in the data6 .

Three other proxies of volatilities may be used: 1) daily realized variance, 2) daily realized volatility in standard deviation form, and 3) daily realized volatility. In what follows, we present the results of measuring the convenience yield in carbon markets using these three proxies.

Here, instead of AR components, we introduce the carbon spot price level and the realized volatility measure for the 2008 carbon futures contract as exogenous regressor of the convenience yield in Eq. ( 3). In Table 5, the specification with carbon spot price levels and realized volatility measures of carbon futures (up to lag 5) provides inferior results for the modelling of the 2008 convenience yield variable compared to AR processes. Indeed, the R 2 is merely equal to 35%.

As shown in Table 6, if we reduce the number of lags to 1 for the realized volatility component, the R 2 falls to 29%.

If we remove the realized volatility component as an exogenous regressor in Eq. ( 3), the specification with carbon spot price levels only has an explanatory power of 25% (Table 7). The latter result reinforces the belief that volatility is an important determinant in the modelling of the 2008 convenience yield variable for carbon spot and futures prices.

Finally, we may use moving averages as another proxy for volatility. In addition to the AR(4) configuration of the convenience yield variable, the inclusion as exogenous regressors of carbon spot price levels and moving averages for the volatility component yield to the best estimates.

In Table 8, the R 2 is indeed equal to 74% for such a regression. The inclusion of realized volatility measures (instead of moving averages) as a proxy for volatility does not yield to superior results (Table 9).

To sum up the results obtained, we find that (i) using an AR(4) process can already explain 74% of the convenience yield for CO 2 spot and futures prices; and (ii) using realized volatility significantly impacts the convenience yield, but the explanatory power is reduced to 35%. The methodology conducted here therefore provides useful information for market players in need to hedge against a potential carbon price risk. Besides, the analysis is based on most recent findings from financial econometrics.

Conclusion

The EU ETS has fostered the development of market place for CO 2 allowances in Europe. The most liquid spot trading place is the BlueNext in Paris, while futures and options trading takes place at the European Climate Exchange in London.

This article models the convenience yield in the European carbon market. The convenience yield stems from differences in spot and futures prices, and can explain why firms hold inventories. Hence, this analysis appears of particular importance for risk management and to traders in Europe.

Compared to previous literature, our approach builds on [START_REF] Borak | Convenience Yields for CO 2 Emission Allowance Futures Contracts[END_REF] by taking explicitly into account various proxies for the volatility of carbon prices that may be of importance for the behaviour of the convenience yield term. More precisely, we focus on daily and realized measures to proxy for volatility of carbon prices, and we assess their respective importance in the modelling of the convenience yield. Highfrequency data for CO 2 allowances futures have been gathered from the European Climate Exchange.

We show that the 2008 convenience yield variable for carbon spot and futures prices may be modelled as a highly persistent variable. Indeed, after testing for various specifications with spot price levels, realized volatility measures or moving averages as exogenous regressors for volatility, the best results are achieved with a simple AR(4) specification.

Besides, we find evidence that daily and realized volatility measures are statistically significant in modelling the convenience yield variable. This result suggests that the relationship between carbon spot and futures prices is not linear.

That spot price levels, moving averages or intraday volatility estimates contribute to explain the convenience yield only at the margin (the explanatory power of regressions with autoregressive components alone is indeed very high) contains some useful information in terms of forecasting for utilities facing the need to hedge against carbon price changes.

First

  Figure 1: Carbon Spot and Futures Prices from April 22, 2005 to January 16, 2009 Source: Bluenext, European Climate Exchange

Figure 2 :

 2 Figure 2: Carbon Spot and Futures Prices (top) and convenience yield (bottom) from January 2 to December 15, 2008 Source: Bluenext, European Climate Exchange, Euribor

Figure 3 :

 3 Figure 3: Volatility signature plot for 2008 carbon futures contract with sampling frequencies ranging from 0 to 80 minutes.

Figure 4 :

 4 Figure 4: Daily realized variance (top panel), daily realized volatility in standard deviation form (middle panel), and daily realized volatility (bottom panel) for the 2008 carbon futures contract

Table 1 :

 1 Descriptive statisics for the daily data used from January 2 to December 15, 2008 Source: ECX, Bluenext, Thomson Financial Datastream Note: Futures denotes the ECX December 2008 carbon futures price, Spot the Bluenext carbon spot price, Euribor the risk-free rate, Spot price level the regressor of the carbon spot price against a constant, RV vol proxy is the measure of realized volatility for the ECX December 2008 carbon futures contract as computed byChevallier and Sévi (2009a), MA vol proxy is another proxy for volatility using moving averages of the carbon price, and Std.Dev is the standard deviation.

		Futures	Spot	Euribor Spot price level RV vol proxy MA vol proxy
	Mean 22.88945 22.65393 4.895124	-0.000724	0.015279	0.010033
	Median 23.63000 23.38000 5.050000	0.000199	0.006450	0.009269
	Maximum 29.33000 28.73000 5.448000	0.023344	0.329379	0.018092
	Minimum 13.72000 13.70000 3.334000	-0.039217	8.17E-06	0.003589
	Std. Dev. 3.591024 3.507532 0.469375	0.010621	0.031766	0.003374
	Skewness -0.810727 -0.805231 -1.392180	-0.659208	6.733913	0.424142
	Kurtosis 2.962105 2.903512 4.612382	3.960755	59.51318	2.138923
	Observations	240	240	240	240	240	240

Table 2 :

 2 ADF Test Statistic for the 2008 convenience yield variable with carbon spot and futures prices

	t-Statistic Prob.*

Table 3

 3

	Variable Coefficient Std. Error t-Statistic	Prob.
	AR(1)	0.368329 0.072681 5.067725 0.0000
	AR(2)	0.189735 0.076136 2.492068 0.0136
	AR(3)	0.214808 0.076070 2.823817 0.0053
	AR(4)	0.190556 0.072244 2.637655 0.0091
	Adj. R 2	0.735475
	AIC -1.179760
	SC -1.110389
	LL	113.7177
	F-Stat.	0.000000

presents estimation results with ARMA filtering. The best specification is obtained with autore-

Table 3 :

 3 AR(4) Test Results for the 2008 convenience yield variable with carbon spot and futures prices Note: The dependent variable is the 2008 convenience yield for carbon spot and futures prices. Std.Error is the standard error, Prob. is the probability value for statistical significance, AR(p) denotes the lags of the autoregressive components. The quality of the regression is verified through the following diagnostic tests: Ad j.R 2 is the Adjusted R-squared, AIC is the Akaike Information Criterion, SC is the Schwartz information criterion, LL is the Log likelihood, and F -Stat. is the p-value of the F-Statistic.

Table 4 :

 4 AR(4) and ARCH Test Results for the 2008 convenience yield variable with carbon spot and futures prices Note: The dependent variable is the 2008 convenience yield for carbon spot and futures prices. Std.Error is the standard error, Prob. is the probability value for statistical significance, AR(p) denotes the lags of the autoregressive components. In the variance equation, RESID(-1) 2 denotes the ARCH term. The quality of the regression is verified through the following diagnostic tests: Ad j.R 2 is the Adjusted R-squared, AIC is the Akaike Information Criterion, SC is the Schwartz information criterion, LL is the Log likelihood, and F -Stat. is the p-value of the F-Statistic.

		Coefficient Std. Error	Prob.
	AR(1)	0.337292 0.092377 0.0003
	AR(2)	0.228832 0.084789 0.0070
	AR(3)	0.210157 0.085976 0.0145
	AR(4)	0.188615 0.073851 0.0106
		Variance Equation
	C	0.014162 0.001197 0.0000
	RESID(-1) 2	0.177775 0.086428 0.0397
	Adj. R 2	0.732061
	AIC	-1.196933
	SC	-1.092877
	LL	117.3148
	F-stat.	0.000000
	Variable Coefficient Std. Error t-Statistic	Prob.
	Constant -0.341998 0.106699 -3.205265 0.0016
	SPOT	0.041766 0.004371 9.555048 0.0000
	LOGSD(-1)	0.012412 0.024008 0.517001 0.6058
	LOGSD(-2)	0.030391 0.025516 1.191031 0.2352
	LOGSD(-3)	0.026523 0.024341 1.089625 0.2773
	LOGSD(-4)	0.047349 0.025377 1.865828 0.0637
	LOGSD(-5)	0.025009 0.024082 1.038525 0.3004
	Adj. R 2	0.348315
	AIC -0.262496
	SC -0.141097
	LL	31.41210
	F-stat.	0.000000

Table 5 :

 5 Realized Volatility with Lag 5 Results for the 2008 convenience yield variable with carbon spot and futures prices Note: The dependent variable is the 2008 convenience yield for carbon spot and futures prices. Std.Error is the standard error, Prob. is the probability value for statistical significance, SPOT is the exogenous regressor for the level of carbon spot prices, LOGSD(p) denotes the lags of the realized volatility measure of the 2008 carbon futures contract in log-transformation. The quality of the regression is verified through the following diagnostic tests: Ad j.R 2 is the Adjusted R-squared, AIC is the Akaike Information Criterion, SC is the Schwartz information criterion, LL is the Log likelihood, and F -Stat. is the p-value of the F-Statistic.

	Variable Coefficient Std. Error t-Statistic	Prob.
	C -0.464655 0.106403 -4.366921 0.0000
	SPOT	0.039174 0.004490 8.725252 0.0000
	LOGSD(-1)	0.071591 0.019751 3.624685 0.0004
	Adj. R 2	0.296707
	AIC -0.207194
	SC -0.155166
	LL	22.26906
	F-stat.	0.000000

Table 6 :

 6 Realized Volatility with Lag 1 Results for the 2008 convenience yield variable with carbon spot and futures prices Note: The dependent variable is the 2008 convenience yield for carbon spot and futures prices. Std.Error is the standard error, Prob. is the probability value for statistical significance, SPOT is the exogenous regressor for the level of carbon spot prices, LOGSD(p) denotes the lags of the realized volatility measure of the 2008 carbon futures contract in log-transformation. The quality of the regression is verified through the following diagnostic tests: Ad j.R 2 is the Adjusted R-squared, AIC is the Akaike Information Criterion, SC is the Schwartz information criterion, LL is the Log likelihood, and F -Stat. is the p-value of the F-Statistic.

Table 7 :

 7 Carbon Spot Price Level Results for the 2008 convenience yield variable with carbon spot and futures prices Note: The dependent variable is the 2008 convenience yield for carbon spot and futures prices. Std.Error is the standard error, Prob. is the probability value for statistical significance, SPOT is the exogenous regressor for the level of carbon spot prices. The quality of the regression is verified through the following diagnostic tests: Ad j.R 2 is the Adjusted R-squared, AIC is the Akaike Information Criterion, SC is the Schwartz information criterion, LL is the Log likelihood, and F -Stat. is the p-value of the F-Statistic.

	Variable Coefficient Std. Error t-Statistic	Prob.
	C -0.565611 0.101884 -5.551507 0.0000
	SPOT	0.035633 0.004434 8.035785 0.0000
	Adj. R 2	0.245863
	AIC -0.188676
	SC -0.155225
	LL	20.49021
	F-stat.	0.000000
	Variable Coefficient Std. Error t-Statistic	Prob.
	SPOT	0.012241 0.005817 2.104200 0.0368
	VOL MA RET -5.071107 8.302313 -0.610806 0.5421
	AR(1)	0.381148 0.074571 5.111182 0.0000
	AR(2)	0.183629 0.078429 2.341334 0.0203
	AR(3)	0.198623 0.078691 2.524105 0.0125
	AR(4)	0.149432 0.074178 2.014496 0.0455
	Adj. R 2	0.741826
	AIC -1.182653
	SC -1.077027
	LL	113.6214
	F-stat.	0.000000

Table 8 :

 8 Moving Averages as Volatility Proxy Results for the 2008 convenience yield variable with carbon spot and futures prices Note: The dependent variable is the 2008 convenience yield for carbon spot and futures prices. Std.Error is the standard error, Prob. is the probability value for statistical significance, SPOT is the exogenous regressor for the level of carbon spot prices, VOL MA RET is the exogenous regressor for a proxy of volatility using moving averages. The quality of the regression is verified through the following diagnostic tests: Ad j.R 2 is the Adjusted R-squared, AIC is the Akaike Information Criterion, SC is the Schwartz information criterion, LL is the Log likelihood, and F -Stat. is the p-value of the F-Statistic.

	Variable Coefficient Std. Error t-Statistic	Prob.
	SPOT	0.009428 0.006072 1.552626 0.1223
	LOGSD -0.011838 0.014946 -0.792022 0.4294
	AR(1)	0.358909 0.073306 4.896050 0.0000
	AR(2)	0.186561 0.076876 2.426772 0.0162
	AR(3)	0.199642 0.076745 2.601353 0.0101
	AR(4)	0.184547 0.073290 2.518040 0.0127
	Adj. R 2	0.736643
	AIC -1.173727
	SC -1.069671
	LL	115.1566
	F-stat.	0.000000

Table 9 :

 9 Counter-Factual Exercise with Realized Volatility Component (instead of Moving Averages) as Volatility Proxy Results for the 2008 convenience yield variable with carbon spot and futures prices Note: The dependent variable is the 2008 convenience yield for carbon spot and futures prices. Std.Error is the standard error, Prob. is the probability value for statistical significance, SPOT is the exogenous regressor for the level of carbon spot prices, LOGSD(p) denotes the lags of the realized volatility measure of the 2008 carbon futures contract in log-transformation. The quality of the regression is verified through the following diagnostic tests: Ad j.R 2 is the Adjusted R-squared, AIC is the Akaike Information Criterion, SC is the Schwartz information criterion, LL is the Log likelihood, and F -Stat. is the p-value of the F-Statistic.

Note differences between implied volatility (extracted from option prices) and realized volatility (computed from intraday data) in the context of convenience yield estimates for the carbon market are left for further research.

Note we do not consider here breakpoint analysis in the time-series of carbon prices. This aspect is left for further research.

Note that, given the interaction between carbon prices and other energy markets[START_REF] Alberola | Price drivers and structural breaks in European carbon prices 2005-2007[END_REF]), we may also include other factors such as oil and gas prices in order to model the convenience yield for carbon prices. This area is left for further research.

This comment applies in the remainder of the article.

Results from the Ljung-Box-Pierce test statistic are not reproduced here to conserve space, but they may be obtained upon request to the authors.

Acknowledgements

I wish to thank Benoît Sévi for fruitful discussions on this topic, as well as ECX for providing the tick-by-tick data. The usual disclaimer applies.


Pauli V Luoma

Cytochrome P450 and gene activationcholesterol elimination and regression of atherosclerosis: author reply to the commentary

published or not. The documents may come

The Editorial commentary by Olavi Pelkonen and Jukka Hakkola "Cytochrome P450, induction and cholesterolwhat are the connections?" [START_REF] Pelkonen | Cytochrome P450, induction and cholesterolwhat are the connections?[END_REF] regarding my review article "Cytochrome P450 and gene activationfrom pharmacology to cholesterol elimination and regression of atherosclerosis" [START_REF] Luoma | Cytochrome P450 and gene activationfrom pharmacology to cholesterol elimination and regression of atherosclerosis[END_REF] published in this journal gives reason for some comments.

The review does not report, as stated in the commentary, that the effects of all the compounds affecting the fate of cholesterol are partially or totally mediated through the action of cytochrome P450 enzymes. Instead, it focuses, as presented, first on the effects of P450s, and second on the effects of gene activation and gene-activating compounds.

The conclusion of the article "P450-enzymes are essential in the physiological maintenance of cholesterol balance" is according to the commentary "correct only at a very general level." There are, however, numerous studies showing that P450 enzymes are specifically involved. P450s are essential in the formation of oxysterols and bile acids, biologically active necessary signaling molecules and metabolites. A defect in CYP27A1 underlies cholesterol accumulation and clinical manifestations such as xanthomas in tendons and brain and premature atherosclerotic disease, and a defect in CYP7A1 leads to cholesterol accumulation and atherogenic hypercholesterolemia.

The commentary states that the simple hypothesis (which is three decades old) linking P450 and induction phenomenon with cholesterol metabolism, transporters and atherosclerosis "still remains a hypothesis." The statement is no longer valid. Cholestyramine, which by specifically inducing hepatic CYP7A1 reduces low-density lipoprotein cholesterol (LDL-C) and, by inducing apolipoprotein (apo) AI synthesis increases apo AI and high-density lipoprotein cholesterol (HDL-C), has been shown to regress coronary atherosclerosis (angiography) [START_REF] Luoma | Cytochrome P450 and gene activationfrom pharmacology to cholesterol elimination and regression of atherosclerosis[END_REF]. Fibrates induce P450s, including ω-hydroxylases, which generate eicosanoid ligands for peroxisome proliferator-activated receptor (PPAR)-α [START_REF] Raucy | Regulation of CYP2E1 by ethanol and palmitic acid and CYP4A11 by clofibrate in primary cultures of human hepatocytes[END_REF][START_REF] Ng | Cytochrome P450 eicosanoids are activators of peroxisome proliferatoractivated receptor α[END_REF]. Fibrate-caused PPAR-α activation upregulates LXR, transporters and cholesterol elimination to apo AI, and an inhibition of the fibrate-caused P450-induction prevents activation of cholesterol-lowering mechanisms [START_REF] Luoma | Cytochrome P450 and gene activationfrom pharmacology to cholesterol elimination and regression of atherosclerosis[END_REF].

Statins induce P450 enzymes, including ω-hydroxylases, in human liver microsomes [START_REF] Luoma | Cytochrome P450 and gene activationfrom pharmacology to cholesterol elimination and regression of atherosclerosis[END_REF][START_REF] Hsu | Regulation of human cytochrome P450 4F2 expression by sterol regulatory elementbinding protein and lovastatin[END_REF]; ω-hydroxylases generate eicosanoids, which are PPAR-α and PPAR-γ activators and also induce the expression of P450 enzymes [START_REF] Ng | Cytochrome P450 eicosanoids are activators of peroxisome proliferatoractivated receptor α[END_REF]. Statincaused PPAR activation upregulates LXR and cholesterol transporters and promotes cholesterol elimination to apo AI and HDL [START_REF] Luoma | Cytochrome P450 and gene activationfrom pharmacology to cholesterol elimination and regression of atherosclerosis[END_REF]. Statin has also been shown to induce-via SREBP-2 activation-the expression of CYP4F2 [START_REF] Hsu | Regulation of human cytochrome P450 4F2 expression by sterol regulatory elementbinding protein and lovastatin[END_REF] and the liver-specific ABCA1 [START_REF] Tamehiro | Sterol regulatory element-binding protein-2 and liver X receptor-driven dual promoter regulation of hepatic ABC transporter A1 gene expression. Mechanism underlying the unique response to cellular cholesterol status[END_REF]. The drug reduced expression of the peripheral type of ABCA1 in vitro in an LXR-dependent manner [START_REF] Tamehiro | Sterol regulatory element-binding protein-2 and liver X receptor-driven dual promoter regulation of hepatic ABC transporter A1 gene expression. Mechanism underlying the unique response to cellular cholesterol status[END_REF] and increased LXR expression both in hepatocytes and leucocytes in vivo [START_REF] Luoma | Cytochrome P450 and gene activationfrom pharmacology to cholesterol elimination and regression of atherosclerosis[END_REF]. LXR activation has been shown to also upregulate the LDL receptor gene, and SREBP-2 activation provides a link between the regulation of CYP4F2 and LDL receptor [START_REF] Hsu | Regulation of human cytochrome P450 4F2 expression by sterol regulatory elementbinding protein and lovastatin[END_REF].

The statement "convincing evidence linking P450s, induction, cholesterol and atherosclerosis is still lacking" does not take into account current knowledge. Deficient P450-hydroxylase activity leads to cholesterol accumulation, and gene-activating agents prevent or regress atherosclerosis They include effective compounds that, via P450 induction, affect the fate of cholesterol and act against the atherosclerotic process. The effects of both endogenous and pharmacologic gene activation have significance: coronary heart disease is the leading cause of death in the world and cerebrovascular disease the second leading cause.

I thank Olavi Pelkonen and Jukka Hakkola for their valuable comments.


J Lars

G Nilsson

Gunnar Lindberg

Gunilla Ringbäck Weitoft

Arne Melander

Gunilla Ringbäck

Aggregated and individual pharmacy record data. Association between antibiotic and antihypertensive drug use

Aggregated and individual pharmacy record data. Association between antibiotic and antihypertensive drug use

In a previous attempt to relate drug use to morbidity we compared "sickness numbers" and dispensed drug volumes in Swedish municipalities and found positive correlations between sickness numbers and dispensed volumes of most drug types. Antibiotics were an exception; higher antibiotics volumes were associated with lower sickness numbers [START_REF] Nilsson | Correlation of drug utilisation and morbidity at the municipality level. High use of antibiotics associated with low use of antihypertensives[END_REF].

Aggregated pharmacy record data showed that municipalities with antibiotics volumes above the national average had lower volumes of cardiovascular (CV) drugs, primarily of antihypertensives [START_REF] Nilsson | Correlation of drug utilisation and morbidity at the municipality level. High use of antibiotics associated with low use of antihypertensives[END_REF], indicating a possible link between infection and hypertension. However, the employed aggregated pharmacy record data gave no information regarding the relation between the dispensed volumes of antibiotics and CV drugs for individuals. The present study aimed at obtaining such information.

The dispensed drug volumes in Sweden (2006) in DDD of antibiotics (ATC code J01) and the CV drugs diuretics (C03), beta-blocking agents (C07), calcium antagonists (C08) and ACE-inhibitors (C09) per 10-year age groups for patients aged ≥ 40 years were extracted from the Swedish Prescribed Drug Register [START_REF] Wettermark | The new Swedish Prescribed Drug Registeropportunities for pharmacoepidemiological research and experience from the first six months[END_REF]. The number of individuals and percentages of the population with anti-biotics or CV drugs were determined. The expected percentage with both antibiotics and CV drugs was calculated as the product of the proportion of the population with antibiotics times the proportion with CV drugs. The ratios of observed percentage/expected percentage were calculated for each age group.

To determine if the use of antibiotics influenced the CV drug volume per individual, we compared for each age group the average volumes of CV drugs (DDDs of C03-C09) for patients with both antibiotics and CV drugs with the average volumes of CV drugs for patients with no antibiotics.

The extracted data comprised 2.35 million individuals or between 347,000 and 597,000 individuals per age group. The ratios of observed/expected percentages of patients with both antibiotics and CV drugs were higher than 1 in four of the five age groups, indicating that more patients than expected had acquired both types of drugs. The average volumes of CV drugs (average DDD/patient) were 5-6% higher per patient treated with antibiotics than for those not treated with antibiotics.

The present study showed that the proportion of the population with both antibiotics and CV drugs was higher than expected and that the groups of individuals treated with antibiotics on average had 5-6% higher volumes of CV drugs than those without antibiotics. Consequently, the use of antibiotics was not associated with lower volumes of CV drugs at the individual level as might be expected from the previous study using aggregated data [START_REF] Nilsson | Correlation of drug utilisation and morbidity at the municipality level. High use of antibiotics associated with low use of antihypertensives[END_REF]. This is in accordance with reported observations where patients with CV disease were likely to be frequent attenders [START_REF] Vedsted | Drug prescription for adult frequent attenders in Danish general practice: a population-based study[END_REF] of health care centres and to receive prescriptions of antibiotics as comedication [START_REF] Tobi | Studying co-medication patterns: the impact of definitions[END_REF].

Although the behaviour of individuals constitutes the characteristics of municipal populations, individual patient drug use data need not be congruent with the drug use data of the municipalities. Factors of importance for individual drug use are likely to be the patient's medical needs and the readiness to consult health care professionals. At the municipal level, factors such as the availability of health care and local therapy traditions may be important.

To summarize, the results of our study illustrate the difficulties in drawing conclusions regarding therapy from drug use data. Aggregated data indicated a possible link between hypertension and infection [START_REF] Nilsson | Correlation of drug utilisation and morbidity at the municipality level. High use of antibiotics associated with low use of antihypertensives[END_REF]. This was not supported by individual-based data. A functional association between antibiotics and CV drugs, implied by the aggregated data, is therefore less likely.

DO00572; No of Pages J. L. G. Nilsson (*) : G. Lindberg : A. Melander NEPI Foundation, Stockholm, Sweden e-mail: lars.nilsson@nepi.net G. Ringbäck Weitoft Epidemiologic Centre, Socialstyrelsen, Stockholm, Sweden


Zemira Cannioto

Giorgio Tamburlini

Federico Marchetti
email: marchetti@burlo.trieste.it

Statins for children? A word of caution

Statins for children? A word of caution

Sir,

Recently the American Academy of Pediatrics (AAP) issued guidelines on screening and treatment of hypercholesterolemia in children [START_REF] Daniels | the Committee on Nutrition (2008) Lipid screening and cardiovascular health in childhood[END_REF]. According to the new indications, statins are now recommended for younger patients (from 8 years instead of 10 as previously recommended) belonging to the same two risk categories addressed in the previous AAP statement: [START_REF] Daniels | the Committee on Nutrition (2008) Lipid screening and cardiovascular health in childhood[END_REF] LDL persistently above 160 mg/dl despite diet therapy, with other risk factors for cardiovascular disease such as overweight, hypertension, smoking, or familial hypercholesterolemia or (2) LDL >190 mg/dl, despite diet therapy, without other risk factors for cardiovascular disease. Based on these indications, a 9-year-old overweight child with LDL, say, at 175 mg/dl after 6 months of appropriate dietetic restrictions, should start treatment with statins.

We think that the evidence to support such an indication is not strong enough and that a word of caution is necessary. Almost all clinical trials of treatment with statins carried out in children have considered patients affected by familial hypercholesterolemia [START_REF] Arambepola | Statin treatment for children and adolescents with heterozygous familial hypercholesterolaemia: a systematic review and meta-analysis[END_REF]. For the relatively small group of patients affected by this condition, a pharmacological intervention is recommended because of the strong correlation with adult-onset cardiovascular disease [START_REF] Wierzbicki | Familial hypercholesterolaemia: summary of NICE guidance[END_REF]. Children with high LDL cholesterol and other risk factors for cardiovascular disease, i.e., the category of patients now included in the risk categories proposed by AAP, represent a far more numerous population group, but in this group of patients, information on long-term efficacy and safety of treatment with statins is scanty. With respect to the efficacy, as children don't have heart attacks, investigators have looked at the effects of statins on endothelial dysfunction and carotid intimal medial thickness, early markers for atherosclerosis in adults. Controlled studies in children show that statins improve these parameters [START_REF] Jongh | Early statin therapy restores endothelial function in children with familial hypercholesterolemia[END_REF][START_REF] Wiegman | Efficacy and safety of statin therapy in children with familial hypercholesterolemia: a randomized controlled trial[END_REF], but there is no evidence of efficacy on prevention of mortality or morbidity related to adult-onset cardiovascular disease. This is a crucial point since the use of statins in primary prevention of cardiovascular disease appears to be controversial in adults, too [START_REF] Moride | Clinical and public health assessment of benefits and risks of statins in primary prevention of coronary events: resolved and unresolved issues[END_REF].

With respect to safety, clinical trials demonstrated a safety profile of statins for a mean time of treatment of 27 weeks (range 6-96) [START_REF] Arambepola | Statin treatment for children and adolescents with heterozygous familial hypercholesterolaemia: a systematic review and meta-analysis[END_REF], which is an extremely short period of time, especially if life-long treatment is proposed, as it is in the AAP statement. The possibility of adverse effects, which may become apparent only after a long exposure to treatment, also considering the greater cumulative dose when treatment starts in childhood rather than in adult life, cannot be neglected.

Long-term treatments also pose the problem of compliance. It is well known that compliance with preventive treatment is much more difficult to obtain than with curative treatment [START_REF] Poluzzi | Adherence to statin therapy and patients' cardiovascular risk: a pharmacoepidemiological study in Italy[END_REF] but surprisingly this issue was not addressed by the AAP statement.

For children belonging to high-risk groups, we believe that further research is needed before we extend the indications for long-term treatment with statins to hypercholesterolemic children and adolescents. Meanwhile, from a public health perspective, we would rather emphasize the need for a more comprehensive approach to address the epidemic of childhood obesity [START_REF] Davis | Recommendations for prevention of childhood obesity[END_REF] and to prevent the related adult-onset risks [START_REF] Baker | Childhood body-mass index and the risk of coronary heart disease in adulthood[END_REF].


D E M Geurts

W J F M Van Der Velden
email: w.vandervelden@hemat.umcn.nl

K M Hebeda

J M M Raemaekers

Richter's syndrome developing in a patient with adult onset Still's disease

Dear Editor, Autoimmune disorders can be a risk factor for the development of lymphoma [START_REF] Smedby | Malignant lymphomas in autoimmunity and inflammation: a review of risks, risk factors, and lymphoma characteristics[END_REF], and at the same time, malignant lymphoma can present with paraneoplastic syndromes consisting of autoimmune phenomena and even full-blown rheumatoid disorders [START_REF] Varoczy | Malignant lymphoma-associated autoimmune diseases-a descriptive epidemiological study[END_REF]. When a patient presents with predominant clinical symptoms of autoimmune disease awareness is necessary in order not to miss an underlying lymphoma. In this paper, we present a case that illustrates the difficulties encountered diagnosing malignant lymphoma when autoimmune phenomena dominate the clinical presentation.

A 53-year old woman was referred to our hospital for a second opinion. Seven months earlier, she had been admitted elsewhere because of intermittent fevers, rigors, night sweats, anorexia, weight loss, myalgia, arthralgia, and a maculopapular rash. Before this, there was no relevant medical history. On physical examination, no pathological lymph nodes, splenomegaly, or hepatomegaly was present. Laboratory investigation revealed a leukocytosis (34.7× 10 9 /l) with 90% neutrophils, elevated C-reactive protein (222 mg/l), liver test abnormalities [alkaline phosphatase (168 U/l), aspartate transaminase (76 U/l), alanine transaminase (49 U/l)], elevated lactate dehydrogenase (910 U/l), and hyperferritinemia (9,409 μmol/l). Additional laboratory tests and serology did not show signs of rheumatoid diseases, sarcoidosis, viral, or bacterial infections (e.g., Epstein-Barr virus (EBV), cytomegalovirus, parvovirus, human immunodefficiency virus, mycoplasma, and borrelia). Computer tomography (CT) revealed two small lymph nodes in the left axilla, but a biopsy failed to acquire representative material. Eventually, the diagnosis of adult onset Still's disease (AOSD) was made because she met almost all criteria for this disease and no signs of a lymphoproliferative disease were found (Table 1) [START_REF] Masson | Comparative study of 6 types of criteria in adult Still's disease[END_REF]. Treatment was started with prednisone (30 mg once daily), resulting in rapid clinical recovery.

August 2006, the patient was referred because lymphocytosis (absolute lymphocyte count, 5.7×10 9 /l) was noticed repeatedly during follow up, not accompanied by physical complaints or signs or cytopenias. After the prednisone was tapered and finally stopped, with immunophenotyping of peripheral blood (CD5, CD19, CD20, CD23, and Ig lambda positive), the diagnosis of a B-cell chronic lymphatic leukemia (B-CLL) could be established, with RAI stage 0, not necessitating therapy. A relation with the prior diagnosis of AOSD was not considered.

In November 2006, she was readmitted because she experienced a recurrence of the same complaints earlier resulting in the diagnosis of AOSD. Again, leukocytosis (26.8×10 9 /l) with predominantly neutrophils (90%) and 8% lymphocytes was seen. A skin biopsy revealed aspecific dermatitis compatible with AOSD. Because of high fever causing headaches, treatment with naproxen and paracetamol was initiated without benefit. Bone marrow biopsy confirmed the earlier diagnosed B-CLL showing small foci of lymphocytic B cells with expression of CD5 and CD23 (Fig. 1a,b). CT now revealed generalized lymphadenopathy and splenomegaly. A surgical excision of an enlarged axillary lymph node was performed (Fig. 1c,d). Prominent paracortical expansion and infiltration of immunoblastic T lymphocytes (Fig. 1e) was seen admixed with large CD20, CD79a, IgA, and partially CD30-positive blastic B lymphocytes (Fig. 1f). The B lymphocytes were negative for CD5, CD23, and EBV. While clonal immunoglobulin gene rearrangements were detected, no clonal rearrangements of the T-cell receptor genes were found, finally confirming the diagnosis of diffuse large B-cell lymphoma (DLBCL). Within the clinical context, this was regarded as a transformation of the B-CLL, the so-called Richter's syndrome, following 10 months after the diagnosis of AOSD. After initiation of combination chemotherapy [cyclophosphamide, doxorubicin, vincristine, and prednisone (CHOP)] and rituximab, the patient recovered rapidly with normalization of the laboratory abnormalities and the lymphadenopathy. Since the completion of eight courses of R-CHOP, she remains in complete remission and without symptoms, now for more than 12 months.

To our knowledge, this is the first case report on Richter's syndrome in a patient with AOSD. Our patient probably had AOSD that was complicated by the develop-ment of a B-CLL, which later transformed to a DLBCL, because at presentation no evidence was found of lymphoma. The presence of a "subclinical" CLL at disease onset cannot be excluded since no immunoflowcytometry studies were performed on the blood at that time. Autoimmune diseases are known to precede the occurrence of malignant lymphoma [START_REF] Smedby | Malignant lymphomas in autoimmunity and inflammation: a review of risks, risk factors, and lymphoma characteristics[END_REF], and there are, although few, reports on AOSD preceding lymphoma [START_REF] Sono | A case of adult-onset Still's disease complicated by non-Hodgkin's lymphoma[END_REF][START_REF] Trotta | Immunoblastic malignant lymphoma in adult onset Still's disease[END_REF]. However, the time interval between the development of the Richter's syndrome and the initial diagnosis of AOSD in this case was rather short, which is 10 months, realizing that the time interval for lymphomas to occur in autoimmune disorders are rather in the range of many years.

On the contrary, CLL and DLBCL can present with autoimmune phenomena, but these consist predominantly of autoimmune cytopenias, such as autoimmune hemolytic anemia, idiopathic thrombocytopenia, and less common, pure red cell aplasia and neutropenia [START_REF] Kipps | Autoantibodies in chronic lymphocytic leukemia and related systemic autoimmune diseases[END_REF]. Polymorphous autoimmune phenomena do occur in lymphoma, and even reports on lymphoma masquerading AOSD exist, although extremely rare, but in most of these reports, T-cell lymphoma seem to be involved [START_REF] Kawasaki | T cell lymphoma masquerading as adult onset Still's disease[END_REF]. An association between B-CLL or DLBCL and polymorphous autoimmune phenomena such as AOSD has not yet been described. Theoretically, it might be possible for B-cell lymphomas to induce polymorphous autoimmune diseases. The role of B lymphocytes extends beyond producing antibodies, as they are antigen-presenting cells and producers of cytokines, activating and altering the function of T lymphocytes, and their role in the pathogenesis of autoimmune disorders mediated by T lymphocyte deregulation has been established [START_REF] Shlomchik | From T to B and back again: positive feedback in systemic autoimmune disease[END_REF]. In our patient, this course of events is however unlikely because symptoms of AOSD clearly preceded the diagnosis of lymphoma.

The diagnosis of AOSD often poses great difficulty because the disease can mimic many other diseases, such as systemic infections, vasculitis syndromes, and malignancies. Therefore, the final diagnosis must not only rely on the presence of several clinical features but also on firm exclusion of other diseases including lymphoproliferative disorders (Table 1) [START_REF] Masson | Comparative study of 6 types of criteria in adult Still's disease[END_REF]. Furthermore, during follow-up, changing, new, or worsening clinical symptoms could result from lymphoma complicating AOSD creating again a diagnostic problem. On pathological examination, enlarged lymph nodes, as seen in 50-60% of AOSD, are difficult to differentiate from lymphoma on pathological examination, especially from T-cell lymphoma [START_REF] Reichert | Lymph node histology simulating T-cell lymphoma in adult-onset Still's disease[END_REF]. In our patient, during follow-up, the occurrence of lymphadenopathy created the same problem. The predominant infiltrate in her lymph node biopsy consisted of a mixture of immunoblastic B and T cells, as reported before in AOSD [START_REF] Quaini | Immunohistological characterization of lymph nodes in two cases of adult onset Still's disease[END_REF], with only a minority of B-cell blasts, first suggesting the differential diagnosis of a reactive lymph node or a T-cell lymphoma. Determining an aberrant B-lymphocyte population by showing clonal immunoglobulin gene rearrangements was essential for the final diagnose of Richter's syndrome.

No strict recommendations exist on the extent of investigations that are required to exclude lymphoma before diagnosing AOSD. Furthermore, currently, no guidelines exist on screening for lymphoma in patients with autoimmune diseases in general, probably because the absolute incidence is low despite a substantial increased risk. In addition, it might be difficult to detect lymphoma in an early stage in patients with multiple aspecific symptoms, as mentioned earlier. Therefore, in practice, to avoid missing underlying lymphomas, extreme vigilance at diagnosis and during follow up is necessary. In case there is a changing or unusual clinical course, especially regarding lymphadenopathy, thorough investigations are obligatory including histopathology and immunophenotyping of blood, bone marrow, and lymph node biopsy. The determination of an aberrant or clonal lymphocyte population is of utmost importance in finally differentiating between autoimmune diseases and lymphoma.

Fig. 1

 1 Fig. 1 Pathological examination of bone marrow and lymph node of the patient. Histology of the CLL in the bone marrow biopsy (a HE ×40, b CD23 ×40) and the transformation in the axillary lymph node (c HE ×20, d HE ×40, e CD3 ×40, f CD20 ×40). The blastic B-cell proliferation (f) is accompanied by immunoblastic T-cells in the background (e)

Table 1

 1 Criteria for diagnosing AOSD

	Criteria
	Major criteria
	Fever of 39°C or higher, lasting one week or longer
	Arthralgia lasting two weeks or longer
	Typical rash*
	Leukocytosis (10.000/mm 2 or greater) including 80 percent of
	granulocytes
	Minor criteria
	Sore throat
	Lymphadenopathy** and/or splenomegaly
	Liver dysfunction
	Negative RF and negative ANA***
	Exclusion criteria
	Infection (especially sepsis and infectious mononucleosis)
	Malignancies (especially malignant lymphoma)
	Other rheumatoid diseases (especially polyarteritis nodosa, and
	rheumatoid vasculitis with extraarticular features)
	Classification of AOSD requires five or more criteria including two or
	more major criteria. Any disease listed under "exclusions" should be
	absent. All criteria are applicable only in absence of other clinical
	explanations.
	a Macular or maculopapular non-pruritic salmon pink eruption usually
	appearing during fever
	b Lymphadenopathy is defined as recent development of significant
	lymph node swelling, and splenomegaly is confirmed on palpation or
	by echogram.

c

Rheuma factor (RF) is serum that must be negative by routine test for the detection of IgM RF, and serum antinuclear antibodies (ANA) must be negative by routine immunofluorescence test.

Open Access This article is distributed under the terms of the Creative Commons Attribution Noncommercial License which permits any noncommercial use, distribution, and reproduction in any medium, provided the original author(s) and source are credited.