Pasquale Laise

Francesca Di Patti

Duccio Fanelli

Marika Masselli

Annarosa Arcangeli

Deterministic and stochastic aspects of VEGF-A production and the cooperative behavior of tumoral cell colony

Keywords: stochastic model, cancer, hypoxia, angiogenesis PACS: 02.50.Ey, 05.40.-a, 82.20.Uv

A model is proposed to study the process of hypoxia-induced angiogenesis in cancer cells. The model accounts for the role played by the vascular endothelial growth factor (VEGF)-A in regulating the oxygen intake. VEGF-A is dynamically controlled by the HIF-1α concentration. If not degraded, HIF-1α can bind to the subunit termed HIF-1β and so experience translocation to the nucleus, to exert its proper transcriptional activity. The delicate balance between these opposing tendencies translates into the emergence of distinct macroscopic behaviors in terms of the associated molecular concentrations that we here trace back to normoxia, hypoxia and death regimes. These aspects are firstly analyzed with reference to the ideal mean-field scenario. Stochastic fluctuations are also briefly discussed and shown to seed a cooperative interaction among cellular units, competing for the same oxygen reservoir.

Introduction

Hypoxia (O 2 tension below 2.5 mmHg) is a hallmark of several types of solid tumor [START_REF] Tatum | Hypoxia: importance in tumor biology, noninvasive measurement by imaging, and value of its measurement in the management of cancer therapy[END_REF]. Growing clinical evidence postulates a correlation between hypoxia and cancer aggressiveness and metastasis [START_REF] Koong | Pancreatic tumors show high levels of hypoxia[END_REF]. Hypoxia may underlie resistance to radiotherapy and can modify the cancer phenotype through gene regulation. A major cellular response to hypoxia is the stabilization of hypoxia inducible factor 1 (HIF-1), a transcription factor that controls cancer progression, and represents a potential target for therapy [START_REF] Semenza | Targeting HIF-1 for cancer therapy[END_REF]. HIF-1 is a heterodimer composed of a oxygen-dependent HIF-1α and a constitutively expressed subunit termed HIF-1β. The transcriptional response to HIF-1 α varies from general effects, such as the up regulation of anaerobic respiration (virtually in all tumors), to the more tissue-specific effects such as angiogenesis [START_REF] Tatum | Hypoxia: importance in tumor biology, noninvasive measurement by imaging, and value of its measurement in the management of cancer therapy[END_REF]. The tumor-driven angiogenesis is therefore a process pivotal in tumor growth and progression, and is mainly related to the HIF-1 dependent secretion of the angiogenic factor Vascular Endotelial Growth Factor -A (VEGF-A) [START_REF] Vaupel | The role of hypoxia-induced factors in tumor progression[END_REF]. Low O 2 tension may shape the cancer phenotype. According to a recently proposed model for cancer progression [START_REF] Fang | Adaptation to hypoxia and acidosis in carcinogenesis and tumor progression[END_REF], the tumor milieu exerts Darwinian selection in favor of cancer cells. The basis for this selection can be two-fold: firstly, the milieu may activate tumor-specific pathways. A second basis for selection may reside in the ability of cells to protect the intracellular compartment from a hypoxic milieu. Given the important role of hypoxia in cancer progression [START_REF] Fang | Adaptation to hypoxia and acidosis in carcinogenesis and tumor progression[END_REF], targeting the cellular responses to hypoxia, including the triggering of VEGF-A secretion and the ensuing angiogenesis, may form an alternative to, or a component of, the current chemotherapeutical practice.

In this paper, a model is proposed to resolve the complex dynamical interplay between the molecular chemical species that participate to the aforementioned process. In particular, we will start by exploring the mean-field scenario stemming from the hypothesized networks of interaction. This formally amounts to operate in the limit for infinite system size, the species concentrations obeying to a closed set of ordinary differential equations. We target our analysis to the asymptotic (equilibrium) regime highlighting the emergence of different attractive (stable) states of the single cell dynamics, which appear to be selectively chosen depending on the chemical parameters involved. We then turn to consider the evolution of three neighbors cells which are effectively coupled via the external environment, and adjust consequently their interior dynamics. The role of fluctuations stemming from the discreteness of the system under scrutiny, is shortly addressed resorting to direct stochastic simulations. Depending on the selected parameters, normoxic condition can develop as an emergent dynamical equilibrium. Allowing for punctual mutations, so to mimic tumor derive, can alter the stability of the system, and eventually result in a competition between distinct cell populations. Indeed, as we shall demonstrate, cooperative mechanisms arise being mediated by the stochastic component of the dynamics. Interestingly, adjacent cells can also profit from the gained ability of an individual cell unit to survive under hypoxic condition, an effect which fades off in the continuum limit.

The paper is organized as follows. Next section is devoted to introducing the key ingredients of the model, here formulated as an ensemble of chemical equations. Then the mean-field version of the model is studied, which in turn corresponds to work in the infinite system size limit. The focus is on a single cell: different dynamical regimes are identified corresponding to distinct choices of the chemical parameters. Stochastic simulations are also performed to test the adequacy of the theory predictions and challenge in silico the role played by finite size corrections. Then in Sec. 5 we turn to discussing a generalization of the proposed model where three cells are made to interact and compete for the oxygen amount. Finally, in Sec. 6, we sum up and draw our conclusions.

A chemical model of hypoxia

Consider now an individual cell and assume the membrane barrier to divide between interior and exterior regions. Selected molecules can migrate through the semi-permeable membrane, being hence transported from the outside to the inside and viceversa. Aiming at elaborating a sound description of the hypoxia cycle [START_REF] Vaupel | The role of hypoxia-induced factors in tumor progression[END_REF], we select a set of candidate molecules which define a close ensemble of interacting elements, according to our schematic representation of the phenomenon under scrutiny. Oxygen molecules (O 2 ) populate the external milieu. They diffuse and eventually reach the membrane where they happen to combine with the hydroxylases (W), which are consequently turned into an active phase, hereafter labeled W a1 . This process is encapsulated into the following chemical equation:

O 2 + W a -→ W a + E 0 ( 1
)
where a is the reaction rate and E 0 stands for the empty case which the oxygen leaves behind when it gets absorbed by the W. The label 0 refers to the chemicals which are populating the outside region, while the molecules confined inside the cell wall are targeted with the subscript i. By introducing the concept of empty spaces (E 0 outside, and E i inside) we impose that the total number of elements (including the empties) is a constant of the dynamics [START_REF] Mckane | Predator-prey cycles from resonant amplification of demographic stochasticity[END_REF][START_REF] Dauxois | Enhanced stochastic oscillations in autocatalytic reactions[END_REF]. Active hydroxylases W a interfere with the HIF-α, yielding to the degradation via the proteasome pathway:

HIF-α + W a b -→ W + E i (2)
Once the HIF-α is being targeted to degradation, the hydroxylases go back to the primitive configuration, waiting for the next oxygen to drive the transition into its active state. The degraded HIF-α molecule is replaced by the empty inner element E i .

Moreover, the HIF-α can occasionally meet the HIF-β subunit. This encounter necessarily takes place within the cell and gives rise to the VEGF i macromolecule2 . The corresponding chemical equation reads:

HIF-α + HIF-β c -→ VEGF i + E i (3)
VEGF i can eventually abandon the cell interior to occupy an external vacancy E 0 , which is hence mutated into a VEGF 0 element. In formulae:

VEGF i + E 0 e -→ VEGF 0 + E i (4)
Furthermore, the VEGF 0 participates to the network of extracellular reactions which underlay the angiogenesis process. In practice, following a complex cascade of nested chemical reactions, the VEGF 0 attracts the oxygen O 2 , which is hence imagined to replace an (external) empty case. This materializes in turn into a simplified vision of a complex dynamic pathway, which has the sole scope of providing a self-consistent entry to the inspected process. The corresponding chemical equation takes the form:

VEGF 0 + E 0 f -→ VEGF 0 + O 2 (5)
Notice that the VEGF 0 molecule is also a product of the reaction, a working hypothesis that we motivate with the empirical observation that VEGF 0 can attract several oxygen units. This process will eventually come to a stop, as follows VEGF degradation which occurs both inside and outside the cells:

VEGF 0 d -→ E 0 (6) VEGF i d -→ E i (7)
Finally, we should also account for the constitutive generation of HIF-α and HIF-β:

E i gα -→ HIF-α (8) E i g β -→ HIF-β (9)
The scenario outlined above is depicted in Fig. 1, where the main molecular actors are also specified3 . It should be noted that the above model conserves the number of molecules which totals in N , a fact that is ultimately related to the presence of the empties species E.

We should emphasize that the proposed model is intimately stochastic, in its original chemical formulation. The inherent stochasticity emerges as an effect of the discreteness of the medium: finite size corrections, also termed demographic noise, result in a endogenous perturbation which can significantly impact the dynamics of the system as compared to the idealized continuum picture, which formally applies in the limit of infinite constituents. To respect the stochastic nature of the problem and simulate it in silico, one can resort to the celebrated Gillespie algorithm [START_REF] Gillespie | A general method for numerically simulating the stochastic time evolution of coupled chemical reactions[END_REF], a method which makes it possible to exactly monitor the time evolution of the discrete population of mutually interacting entities. The Figure 1: Cartoon for a single cell dynamics as follows the proposed scheme. Oxygen molecules (circles laying outside the cell) are captured by the hydroxylases sitting on the membrane, and so trigger the HIF-α (dark circles, red online, stored inside the membrane) elimination pathways (dashed circle). At the same time, HIF-α can combine with HIF-β (light circles, orange online, stored inside the membrane), yielding to the (barrel like, in the picture) VEGF i . This latter can exit the cell wall, so ensuing in VEGF 0 and drive the incoming oxygen flux. scheme captures in fact the probabilistic nature of the microscopic couplings and enables one to resolve the contribution associated to finite size fluctuations. As we shall demonstrate, stochastic fluctuations can eventually materialize in complex dynamics, which call for sound biological interpretation.

As anticipated, and opposed to the probabilistic approach, one can invoke the deterministic approximation, which in turn corresponds to assuming an infinite population amount and consequently disregarding finite size correlations. This defines the mean-field, continuum level of description, to which next sections are entirely devoted. More specifically, we will start by considering the case of a single cell, which we will discuss with reference to its mean field analogue. Notice that the deterministic formulation results from a straightforward derivation, that moves from the underlying master equation for the stochastic process, as defined by the above set of chemical equations. In the following, we will allude to such procedure, recalling its general philosophy, and just provide the end result for the case under inspection, without insisting on the technicalities. Working within the continuum representation we will then elaborate on the asymptotic dynamics of the recovered system of ordinary differential equations, by identifying distinct fixed points and discussing their associated stability characteristics. The inherent stochastic effects stemming from finite N corrections, will be relegated to section 5.

Master equation

At time t, the status of the chemical system, as ruled by the above reactions, is completely specified once the species amount are assigned. These are discrete entries, one for each species belonging to the ensemble of interacting population, which define the scalar components of the status vector n. Recall that the system is intrinsically stochastic, its evolution being a suite of events which realizes with a prescribed probability. In mathematical terms, this fact implies dealing with the probability that the system is seen in the state n at time t, hereafter denoted P (n, t). This latter quantity adjusts in time and obeys to the master equation:

d dt P (n, t) = n [T (n|n )P (n , t) -T (n |n)P (n, t)] (10)
where T (n |n) denotes the transition rates between two adjacent states, from n to n , compatible with the chemical constraints imposed by equations ( 1)-( 9). For the sake of compactness, we omit here to list the specific forms of the transition rates relative to the model under analysis.

The partial differential equation ( 10) is an exact formulation, which can be in principle analyzed to track the system dynamics. In practice however, it is hard to analytically handle and further progress has to rely on direct numerical investigations.

The Gillespie algorithm represents in particular a viable strategy to simulate the time evolution of the probability distribution P (n, t). The method sets up a Monte Carlo procedure to determine the next reaction that is selected to occur, among all possible ones, and the time τ , when the event takes place. Both selections reflect the assigned transition probabilities, which scales proportional to the number of substrate molecules times the associated reaction constants. Then, based on the selected reaction, the update of the molecule count follows, and time t is replaced by its value increased by the stochastic time interval τ , namely t → t + τ .

Alternatively, the master equation ( 10) can be studied via perturbative methods, aimed at extracting both the idealized continuum picture and the successive finite size corrections. The van Kampen large N expansion is in this respect a widely adopted technique, extensively characterized in the literature [START_REF] Gardiner | Handbook of Stochastic Methods[END_REF][START_REF] Van Kampen | Stochastic preocesses in Physics and Chemistry[END_REF]. In the following section, we shall operate within the van Kampen framework to derive a closed system of ordinary differential equations for the coupled evolution of the concentrations of chemical system (1)-( 9) in the limit of large (formally infinite) N .

Mean-field solution of the single cell model

Let us focus on a single cell, as pictorially represented in Fig. 1. Consider for instance the O 2 species. According to van Kampen [START_REF] Van Kampen | Stochastic preocesses in Physics and Chemistry[END_REF], one can put forward the following ansatz for the associated mean-field (continuous) concentration here labeled O:

O = n O2 N + ξ √ N
where the extensive quantity n O2 is discrete: it is one of the components of the vector n and stands for the number of oxygen molecules. ξ is instead a stochastic variable and contributes with an additive terms which scales as 1/ √ N , so materializing in a perturbation to the idealized mean-field dynamics. Stochastic effects can be hence safely neglected when performing the continuum limit [START_REF] Gardiner | Handbook of Stochastic Methods[END_REF][START_REF] Van Kampen | Stochastic preocesses in Physics and Chemistry[END_REF], the relic concentration obeying to a deterministic set of coupled differential equations. We recall however that beyond the idealized scenario, finite size corrections do matter when real systems are concerned, the intimate discreteness being a crucial ingredient of the dynamics. We shall return on this issue in Section 5, to highlight the important impact of such an endogenous source of stochastic perturbation with reference to the scrutinized model.

For simplicity, by extending the above reasoning to the other involved species, we will here denote by W, W a , V 0 , V i , E 0 , E i , H α and H β the average (meanfield) concentration of W, W a , VEGF 0 , VEGF i , E 0 , E i , HIF-α and HIF-β respectively. According to this notation, the mean-field equations take the form4 :

d dt O = -aOW + f V 0 E 0 (11)
d dt W = -aOW + bW a H α (12)
d dt W a = aOW -bW a H α (13)
d dt H α = -bW a H α -cH α H β + g α E i (14) d dt H β = -cH α H β + g β E i (15) d dt V i = cH α H β -eV i E 0 -dV i (16) d dt V 0 = eV i E 0 -dV 0 (17) d dt E i = bW a H α + cH α H β + eE 0 V i + dV i -(g α + g β )E i (18)
with the conservation law

O + W + V 0 + V i + E 0 + W a + E i + H α + H β = 1 (19)
Notice that the same result could be recovered without the use of the van Kampen's expansion. Indeed, multiplying both sides of the master equation by a scalar component of the status vector n, summing over all n, shifting some of the involved sums by ±1 and neglecting the correlations (which is legitimate in the N infinite limit) one recovers the same mean field system.

From Eq. ( 12) and ( 13), it straightforwardly follows that W a + W stays constant: This invariant quantity is hereafter set equal to n 0 and specified by the initial condition. To investigate the dynamics of the above system ( 11) -(18), and in particular elaborate on its asymptotic evolution, one can look for stationary solutions. These are the fixed points of the dynamics, which can be explicitly found by setting to zero the derivatives in ( 11) -( 18) and solving the obtained algebraic system, where the (final) concentration amounts enter as unknowns.

Stationary states

Solving for the fixed points returns three possible solutions. These latter are here characterized and labeled via the progressive indexes i = 1, 2, 3, associated to the predicted asymptotic concentrations. Two solutions can be readily obtained via direct inspection of the above equations, and correspond to opposite tendencies that we will hereafter call normoxia and death respectively.

Normoxia conditions are assumed to correspond to the asymptotic solution:

V 0,1 = V i,1 = H α,1 = W 1 = E i,1 = 0. (20)
HIF-α is turned off and so are the inner/outer concentrations of VEGF. The HIF-α degradation mechanism hence prevails over the alternative pathway that yields VEGF production, and which could in principle contribute to enhance the oxygen supply. The amount of O 2 , here quantified by O 1 , is on the contrary sufficient to guarantee the correct cell functioning. Hydroxylases are frozen in their active state W a and consequently stimulate the repression of the HIF-α quota 5 . From the conservation relations we immediately deduce the following relation:

W a,1 = n 0 (21) O 1 + E 0,1 + H β,1 = 1 -n 0 (22)
The outer cases are shared between oxygen and empties, as exemplified by the latter equation, in a proportion that cannot be predicted a priori.

As opposed to normoxia condition, a trivial solution of the mean field problem is found when:

W a,3 = V 0,3 = V i,3 = H α,3 = E i,3 = O 3 = 0 (23)
with in addition W 3 = n 0 and E 0,3 + H β,3 = 1n 0 . In this case, the cell is dead: the oxygen amount is zero and the hydroxylases cycle is not active, so denouncing the lethargic state of the cell device. An intermediate asymptotic stationary state is also possible, when the cell balances the oxygen consumption by properly adjusting the VEGF production cycle, at equilibrium. Upon manipulation of the fixed point equations associated to the system ( 11) -( 18), the equilibrium solution for E 0,2 is found to verify the following quadratic equation:

f e d E 2 0,2 - g α g β -1 dE 0,2 -d g α g β -1 = 0 (24)
The previous equation admits a solution only if g α > g β , otherwise is reduced to a sum of three positive terms which cannot mutually cancel for positive coefficients and state variables. When g α = g β , then E 0,2 = 0, which corresponds to having the largest allowed oxygen concentration O 2 = 1n 0 . This latter can alternatively be seen as a limiting normoxia condition. All externally available cases are simultaneously occupied by oxygen molecules, which can still effectively contrast the HIF-α generation, while providing the necessary supply for the cell functioning (W a,1 = 0). Beyond this condition, a quota of HIF-α survives the degradation and sustains the VEGF production process. VEGF molecules transported outside the cell support the afflux of oxygen and so contribute to stabilize the quota of HIF-α. In practical terms, the condition g α = g β (when the HIF-α and HIF-β production rates equalize) signals a transition between two distinct asymptotic regimes, a fact that could be further appreciated by quantifying the values assumed by the remaining involved concentrations. We identify the intermediate fixed point with the hypoxia conditions: The cell is in fact alive thanks to a dedicated feedback that integrates the original oxygen amount via the VEGF production. Once g α > g β the solution of Eq. (24) reads:

E 0,2 = d e(g α -g β ) + e(g α -g β ) [e(g α -g β ) + 4f g β )] 2f eg β (25)
The equilibrium values of the other variables cannot be univocally determined. However, given the value of E 0,2 as estimated above, the following ratios are fixed:

O 2 W 2 V 0,2 = f a E 0,2 W a,2 H α,2 V 0,2 = f b E 0,2 H α,2 H β,2 V 0,2 = g β c f g α -g β E 0,2 E i,2 V 0,2 = f g α -g β E 0,2 V i,2 V 0,2 = d e 1 E 0,2
Notice that an implicit relation between W a,2 and V 0,2 is also derived by making use of Eq. ( 19). More interestingly, for the case g α > g β , E 0,2 can be shown to increase with g α , while keeping all the other parameters constant. The maximum value that can be eventually attained is E 0,max = 1n 0 , a limiting condition that, together with g α > g β , identifies the portion of parameters plane for which the hypoxia fixed point do exist. More specifically, the largest the values of f , i.e. the rate of oxygen production via the VEGF cycle, the wider the window in g α that delimits the hypoxia regime. Working in the (g α , f) plane we can analytically determine the boundaries of the domain of existence of the non trivial fixed point. Its equation can be deduced by imposing E 0,2 = E 0,max into Eq. (24) which immediately yields to

f = edE 0,max + d 2 g β eE 2 0,max g α - edE 0,max + d 2 eE 2 0,max
This is a line which crosses the axis f = 0 when g α = g β , see Fig. 2 where three domains, respectively I, II, III, are outlined. Notice that the normoxia and death solutions can in principle extend over all the parameters plane, while the hypoxia is solely bound to region II. The hypoxia region shrinks as f is reduced and eventually fades off when the VEGF's ability to carry oxygen is hindered. We recall in fact that within our simplified formulation, the oxygen income is ultimately controlled by the parameter f , and consequently linked to the available VEGF 0 quota. It is hence plausible that that region here identified to correspond to the hypoxia regime vanishes when f is set to zero. Numerical simulations indicate that the basin of attraction of the normoxia state prevails in region I, while the death configuration dynamically takes over in region III.

In region II, it is the hypoxia condition to represent the stable attractor of the dynamics. For this reason, region I, II and III are respectively termed normoxia, hypoxia and death, and the caption of Fig. 2 reflects such a choice. To gain further insight into this empirical observation, we turn in the next section to studying the stability properties of the identified fixed point. The analysis is carried out via combined numerical and analytical means, as the complexity of the model prevents us from performing a complete analytical treatment. The existence of the normoxia (I), hypoxia (II) and death (III) regions is displayed in the (gα, f) plane, as predicted by the mean-field calculations for N = 1000, n 0 = 20, e = 0.9, d = 0.4 and g β = 0.5. The transition from I to II is found to occur when gα = g β , i.e. when the HIF-α production rate prevails over the HIF-β one. Then the exceeding HIF-α is stabilized to an asymptotic fixed concentration via the additional VEGF pathway. It can be further appreciated that region II gets compressed when reducing the control parameter f .

Stability analysis

In this section we set down to study the stability properties of the three families of equilibrium points, as identified above. To this end we calculate the Jacobian matrix associated to the system (11)-(18), explicitly given in the Appendix and evaluate the eigenvalues relative to the three fixed points, namely the normoxia, the hypoxia and the death.

Consider first the solution that corresponds to the cellular death, previously labeled with the index 3 in the concentration. We recall that W a,3 = V 0,3 = V i,3 = H α,3 = E i,3 = O 3 = 0, W 3 = n 0 , while the concentration E 0,3 and H β,3 can in principle take any value, constrained by the relation E 0,3 + H β,3 = 1-n 0 , which ultimately specifies the family of possible solutions. The associated characteristic polynomial reads:

λ 3 (λ + an 0 )(λ + d) det ⎛ ⎝ λ + cH β,3 0 -g α -cH β,3 λ + eE 0,3 + d 0 -cH β,3 -eE 0,3 -d λ+ g α + g β ⎞ ⎠ = 0 (26)
where λ stands for the eigenvalue. It is evident that λ 1 = λ 2 = λ 3 = 0, λ 4 = -an 0 and λ 5 = -d are possible solution of Eq. (26). The remaining three eigenvalues cannot be analytically determined. For g α > g β (a condition that embraces region II and III of the parameters plane depicted in Fig. 2), analytical progress is possible: just one stable solution exists among those of type 3 that are in principle admissible. To substantiate this observation, let us denote with J sub 3 the non trivial sub-matrix of the full Jacobian (26):

J sub 3 = ⎛ ⎝ -cH β,3 0 + g α +cH β,3 -eE 0,3 -d 0 +cH β,3 +eE 0,3 + d -g α -g β ⎞ ⎠ .
It is easy to see that det(J sub 3 ) = cH β (g αg β )(d + eE 0,3 ), and so det(J sub 3 ) > 0 ∀H β,3 > 0 and g α > g β . Positive determinant implies unstable solutions, as can be immediately appreciated with simple algebraic considerations. This means that the fixed point corresponding to the death state is always unstable when g α > g β , and H β,3 = 0, this latter quantity being positively defined. At variance, stable solution can manifest in the region specified by condition g α > g β , when H β,3 = 0. This latter request equivalently implies E 0,3 = 1n 0 , and removes the degeneracy (at least in region I and II) by uniquely characterizing the fixed point of type III selected by the dynamics. The only stable fixed point of type "death" in region II and III, is now completely specified and the three lacking eigenvalues can be readily calculated, resulting in λ 6 = 0, λ 7 = -e(1n 0 )d and λ 8 = -g αg β . Solutions of type 3 are clearly stable also in region I, suggesting that the cell can always meet the conditions that drive its death, an asymptotic fate which is ultimately determined by the selected initial datum (corresponding to the environmental conditions).

However, direct numerical inspection, suggests that within region II, the system tends to converge to the hypoxia solution, this latter being hence dynamically favored over the stable death regime. Tracing the respective basins of attraction can be tackled via a dedicated campaign of numerical simulations. Biologically, it can be hence speculated that for sufficiently large values of the parameter g α the cell opposes the death derive, by compressing its respective basin of attraction upon activation of the hypoxia cycle.

In region I, when g α < g β , an infinite family of solutions of type 3 can in principle manifest, depending on the specific initial datum. The associated eigenvalues can be calculated numerically by tuning both H β,3 and E 0,3 within the allowed range of variation, pointing to the stability of the solutions. Again, the cell can always die in response to particularly stressing conditions. To complete the description of region I, one needs to address the stability properties of the normoxia solution. In this case, repeating the above strategy, the characteristic polynomial turns out to be

λ 3 (λ+aO 1 )(λ+d) det ⎛ ⎝ λ + bW a,1 + cH β,1 0 -g α -cH β,1 λ + eE 0,1 + d 0 -bW a,1 -cH β,1 -eE 0,1 -d λ+ g α + g β ⎞ ⎠ = 0
which immediately gives λ 1 = λ 2 = λ 3 = 0, λ 4 = -aO 1 and λ 5 = -d. In this case the determinant of the residual 3 × 3 sub-matrix does not return any sensible information on the sign of the real part of the three unknown eigenvalues.

To overcome this limitation we proceed as follows: we first integrate numerically the governing system of differential equations and then select the attained equilibrium fixed point, as an entry to calculate the full set of eigenvalues. Among the missing eigenvalues, two (which can be complex or real depending on g α ) always display a negative real part, while the third one is real and can change the sign as g α is increased. This situation is illustrated in the inset of Fig. 4: the real eigenvalue is negative over a range of parameters that practically corresponds to region I, and it approaches zero for g α g β . This result holds true for different parameters choice (i.e. the value of f ) and selected initial conditions. It can hence be argued that the normoxia condition is stable within region I and that it gets unstable when the boundary of region II, the domain of existence of the hypoxia solution, is (continuously) approached. Indeed, within region I, the competition between the two attractors of type 1 (normoxia) and 3 (death) is dynamically resolved to favor the former, as confirmed by direct numerical inspection.

In conclusion, the proposed model can alternatively yield to three distinct dynamical regimes, depending on the specific choice of the chemical parameters. The normoxia condition is stable in region I and loses its stability when the edges of the domain II are approached. In this latter region it is the hypoxia solution to take over, by prevailing over the death regime which instead dominates in region III. Interestingly, death is always an attractive state of the dynamics, irrespectively of the selected parameters, a tendency that the cell efficiently opposes in region I and II. Formally, it is hence not accurate to talk about a transition among the three states, the landscape of possible stable attractors being more complex as outlined above. However, to keep the message simple, we shall hereafter reflect the practical segregation into three distinct zones, as virtual transitions. Notice that a punctual mutation of the reference parameters could substantially alter the cell functioning and determine its death. In the next subsection we turn to numerical simulations to validate the proposed scenario. The simulations here performed refer both to the idealized mean-field dynamics and to the stochastic scenario which respects the inherent discreteness of the chemical medium.

Numerical simulations vs. mean field theory

In this section we report about a campaign of simulations aimed at validating the theoretical picture depicted above. Mean-field ordinary differential equations ( 11)-( 18) can be straightforwardly integrated via a standard Runge-Kutta scheme. In Fig. 3 the time evolution of VEGF i and active hydroxylases is reported (solid lines) for different values of the parameters. Panels (a) and (b) refer to a choice of the parameters that yield to the normoxia condition: Once the transient has died out, the system achieves its asymptotic state where VEGF molecules are absent. At variance, W a attains its equilibrium concentration, different from zero. The same applies to oxygen (not displayed here). When increasing the parameter g α , which controls the HIF-α production rate, a transition occurs towards what we termed the hypoxia regime, see panels (c) and (d). Here the HIF-α quota is sensibly different from zero, implying the existence of a positive feedback reaction on the oxygen production, via the VEGF molecules. By further augmenting g α , one gradually approaches the boundary of the third allowed region, as previously outlined: When g α is large enough, all chemical concentrations drop to zero, including the W a amount, which is here assumed to signal the cell activity. This scenario is clearly displayed in panels (e) and (f) of Fig. 3.

Beside integrating the system of differential equations ( 11)-( 18), one can turn to direct stochastic simulations, which enables one to treat the interacting molecules as diffusing massive entities, so respecting the finiteness of the chemical environment. To this end, we resort to the Gillespie algorithm: as already recalled, this is a sophisticated Monte Carlo scheme which enables to inspect the actual dynamics of an ensemble made of discrete actors, whose relative interactions are specified by an assigned set of chemical rules. The outcome of the stochastic integration is superposed in Fig. 3, to the curves determined within the mean-field approximation. The intrinsic stochastic noise materializes as a local disturbance of the reference mean-field profile, suggesting that the discreteness of the environment contributes modestly to the system dynamics. As opposed to this vision, we will prove in the following that finite size correction can be important and can eventually result in drastic modification of the dynamical behavior, as predicted within the mean-field scenario.

Before turning to analyze these aspects, we provide in Fig. 4 a global picture of the distinct dynamical regimes that can be faced when tuning the (f, g α ) parameters: Here the asymptotic concentration of VEGF i , as calculated upon integration of Eqs. ( 11) -( 18) is plotted with a colorcode. The existence of three different regions, is clearly confirmed, the transition from normoxia (I) to hypoxia (II) being sharper and hence more evidently appreciated. As a further complement, we enclose a section of the scanned parameters plane, by reporting the recorded concentration of VEGF i versus g α , at fixed f . The double transition is shown, bearing intriguing analogy with the phenomenon of re-entrant phase transitions.

Dynamics in a multicellular environment: Synchronization and cooperative effects

To shed light onto the cooperative effects that might eventually arise in a colony of interacting cell units, we here consider N c = 3 cells obeying to the chemical equations hypothesized above. Now the oxygen fills a shared reservoir, from which any individual cell can draw. The cells are hence mutually, though indirectly, communicating via the oxygen amount. A schematic layout of the interaction network is depicted in Fig. 5. In this case we need to introduce an extra label to specify the cell to which the molecules belong to. In practical terms, one is forced to deal with the molecules VEGF j i , E j i , HIF-α j and HIF-β j with j = 1, . . . , N c . The governing chemical equations follow as a straightforward generalization of Eqs. (1)- [START_REF] Mckane | Predator-prey cycles from resonant amplification of demographic stochasticity[END_REF]. The mean-field system calculated from the underlying master equation is now composed by 4 × N c + 3 equations (not explicitly given here), the coupling being indirectly provided by the species hosted in the outer environment which constitutes a shared reservoir.

Working within this generalized setting, interesting cooperative effects arise, as supported by the discrete component of the dynamics. In the finite N regime, correlations exist between the evolving species and translate into important contributions that can affect the system dynamics. The correlations are instead neglected within the mean-field approximation, formally valid when the container is filled with an infinite amount of molecular constituents. To demonstrate the importance of finite size corrections, we assign to two cells a set of parameters (including the g α value) which is found to yield to the death condition, in the mean field regime and according to the preceding classification. The third cell is instead imagined to bear a g α parameter in the hypoxia domain (while keeping the other parameters unchanged), as identified above. What is going to happen to the considered system of three cells?

In Fig. 6 the evolution of hydroxylase is reported as a function of time. The solid lines (red online) stand for the mean-field dynamics, obtained by numerical integration of the generalized set of ordinary differential equations which accounts for the three inspected cells. Two cells are destined to die (no residual W a ), while the third survives (W a are eventually present). Interestingly the stochastic dynamics is peculiarly different, as displayed in Fig. 6 (wiggling curves, cyan online). As resulting from a cooperative mechanism, ultimately driven by the finiteness of the simulated medium, and so not captured within the corresponding continuous description, the two cells destined to death are rescued by the third cell. This latter can sustain the angiogenesis process and so integrate the oxygen quota which is necessary for the correct functioning of its nearest neighbors. The rescued cells do partially contribute by stimulating the production of VEGF. Synchronized oscillations are observed in the recorded dynamics of the active hydroxylases, as displayed in the zoomed inset of Fig. 6.

In conclusion, and based on the above findings, it is tempting to speculate that the observed degree of cooperation turns out particularly crucial for understanding the adaptation of tumoral cells in a changing environment. Under stressing condition, it is well known that tumoral cells can undergo mutation, changing their specific abilities and fully exploiting their inherent flexibility. Assume that under stressing condition the oxygen income is not enough to guarantee the cell functioning. Then the tumoral cells would undergo a critical phase, which could induce punctual mutations. In particular, imagine that following such random adjustments, a cell enters the regime of hypoxia (in our model, by properly tuning the corresponding g α parameter). This corresponds to reducing the HIF-α production rate, so consequently freeing a quota of oxygen, otherwise employed in the HIF-α degradation pathway. Therefore not only the tumoral cell interested by the mutation can survive, but also the cells populating the local neighborhood (which are plausibly also belonging to the tumoral family although not necessarily affected by the same mutation), follows a cooperative interaction. This effect is here reproduced within a sound dynamical picture where finite size effects prove to be crucial.

Conclusions

In this paper we have discussed a simplified model to quantitatively address the dynamical mechanisms which ultimately underlie the angiogenesis process.

The dynamical interplay between a limited number of molecular species is resolved and investigated both via numerical and analytical techniques. The level of complexity of the proposed description results from a trade-off balance between the request of incorporating a fair amount of biological information, and the need of making the problem analytically tractable. Within the proposed formulation, the cell can follow different routes, depending on the specific values of the involved chemical parameters. Distinct macroscopic regimes are here discussed and associated to normoxia, hypoxia and death conditions respectively. The rate of HIF-α production, as well as the parameter that controls the ability of VEGF to call for new oxygen, are shown to play a crucial role in determining the transition between the admissible dynamical regimes. These latter are analytically explained via a mean-field based calculation, which is fully confirmed by direct numerical simulations. The role of stochastic noise is also elucidated and shown to drive a degree of synchronization between the cells belonging to a given colony. This is a spontaneously emerging feature, stemming from the intimate discrete nature of the molecular environment, which can potentially relate to the attested adaptivity of localized families of tumoral cells. In fact, it is in general believed that cancer cells can react to the external stressing stimuli (as for the lack of oxygen) by favoring a specific mutation, which renders them more adaptable to the hosting environment. While this is a customarily evoked scenario, the simultaneous occurrence of the same mutation in all the cancer cell population appears a highly improbable event. Cooperative survival mechanisms, as those here identified, could in turn contribute to explain the observed robustness, via a self-consistent dynamical argument. As an additional point, the results of this paper can be of interest in studying the growth of a tumor cord characterized by an aerobic-anaerobic switch triggered by hypoxia [START_REF] Astanin | Mathematical modelling of the warburg effect in tumour cords[END_REF].

Acknowledgments

The work was supported by grants from the Associazione Genitori contro le Leucemie e Tumori Infantili Noi per Voi, Associazione Italiana per la Ricerca sul Cancro, Istituto Toscano Tumori to AA, Ente Cassa di Risparmio di Firenze to the Dipartimento di Patologia e Oncologia Sperimentali, Università degli Studi di Firenze.

Appendix A. Jacobian matrix

The Jacobian matrix associated to the system (11)-(18) reads VEGF i and Wa) are depicted for gα = 0.0017, while keeping the other parameters unchanged. A residual quota of VEGF i molecules is present, which sustains the hypoxia pathway as described in the main text. Finally, when further increasing gα (gα = 0.017 in panels (e) and (f), respectively referring to the VEGF i and Wa concentrations), the cell dies, no activity being asymptotically detected. Solid lines refer to the integration of the mean-field dynamics ( 11) -(18), while the associated irregular profiles stand for the corresponding stochastic simulations. For all the plots N = 10000, the initial number of O 2 , W, W * and HIF-α is 1000, that of HIF-β, VEGF i , VEGF 0 and E i is equal to 100 and, hence, the initial number of E 0 is 5600. The data follows from a direct integration of the mean-field Eqs. ( 11) -( 18) and are represented using the colorcode displayed in the legend. Right panel: The concentration of VEGF i is reported as a function of gα, for f = 0.005 (the other parameters are set as in Fig. 3. A double transition is clearly displayed, the VEGF i acting as a stabilizer of the cell dynamics for a compact range of gα values. Inset a non trivial (see text) real eigenvalue associated to the first family of equilibrium points (normoxia), as a function of gα. The eigenvalue is seen to approach zero, i.e. the threshold of stability point, when gα = g β . The eigenvalue is calculated via an ad hoc semi-analytical procedure which is described in the main body of the paper. Figure 6: Wa evolution as a function of time, for a system of three cells, sharing the same oxygen reservoir. The parameters are set so to have two cells in the mean-field region deputed to death, while the third cell is made to function in the hypoxia domain, via a appropriate mutation of the parameter gα. All other parameters are identical for the three cells. Solid lines (red online) represent the mean-field dynamics, while the irregular curves (cyan online) stand for the stochastic simulations. Inset: a zoom of the dynamical evolution is represented.

J = ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ -aW -aO 0 0 0 0 f E 0 0 -aW -aO +bH α +bW a 0 0 0 0 +aW +aO -bH α -bW a 0 0 0 0 0 0 -bH α -bW a -cH β -cH α 0 0 g α 0 0 0 -cH β -cH α 0 0 g β 0 0 0 +cH β +cH α -eE 0 -d 0 0 0 0 0 0 0 +eE 0 -d 0 0 0 bH α bW a + cH β cH α eE 0 + d 0 -g α -g β ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ 0

Figure 2 :

 2 Figure2: The existence of the normoxia (I), hypoxia (II) and death (III) regions is displayed in the (gα, f) plane, as predicted by the mean-field calculations for N = 1000, n 0 = 20, e = 0.9, d = 0.4 and g β = 0.5. The transition from I to II is found to occur when gα = g β , i.e. when the HIF-α production rate prevails over the HIF-β one. Then the exceeding HIF-α is stabilized to an asymptotic fixed concentration via the additional VEGF pathway. It can be further appreciated that region II gets compressed when reducing the control parameter f .

Figure 3 :

 3 Figure 3: Simulated dynamics of system (11) -(18) for different choices of the chemical parameters. Panels (a) and (b) report the time evolution of the concentration of VEGF i and Wa respectively, for a=0.01, b=0.01, c=0.05, e=0.01, f =0.0005, d=0.0001, gα=0.000002, g β =0.0008. This is the normoxia condition, when the cell functions, beyond the transient, without activating the hypoxia cycle. In panel (c) and (d) the homologous quantities (resp. VEGF i and Wa) are depicted for gα = 0.0017, while keeping the other parameters unchanged. A residual quota of VEGF i molecules is present, which sustains the hypoxia pathway as described in the main text. Finally, when further increasing gα (gα = 0.017 in panels (e) and (f), respectively referring to the VEGF i and Wa concentrations), the cell dies, no activity being asymptotically detected. Solid lines refer to the integration of the mean-field dynamics (11) -(18), while the associated irregular profiles stand for the corresponding stochastic simulations. For all the plots N = 10000, the initial number of O 2 , W, W * and HIF-α is 1000, that of HIF-β, VEGF i , VEGF 0 and E i is equal to 100 and, hence, the initial number of E 0 is 5600.

Figure 4 :

 4 Figure4: Left panel: The asymptotic concentration of VEGF i is plotted for different values of (gα, f). The data follows from a direct integration of the mean-field Eqs. (11) -(18) and are represented using the colorcode displayed in the legend. Right panel: The concentration of VEGF i is reported as a function of gα, for f = 0.005 (the other parameters are set as in Fig.3. A double transition is clearly displayed, the VEGF i acting as a stabilizer of the cell dynamics for a compact range of gα values. Inset a non trivial (see text) real eigenvalue associated to the first family of equilibrium points (normoxia), as a function of gα. The eigenvalue is seen to approach zero, i.e. the threshold of stability point, when gα = g β . The eigenvalue is calculated via an ad hoc semi-analytical procedure which is described in the main body of the paper.

Figure 5 :

 5 Figure5: Cartoon of the interaction scheme of three cells. Oxygen molecules (circles laying outside the cells' walls) are now shared, and provides an indirect couplings among the three cells. Each cell is then supporting a whole cycle of internal reactions, as already exemplified in Fig.1. For the symbols legend, see caption of Fig.1

Hydroxylases are intracellular enzymes, organized in three families respectively PHD1, PHD2, PHD3. Hydroxylases hence populate the internal flow and eventually combine with the oxygen that have penetrated the cell. As formalized in the following, we shall here focus on a more simple picture, and imagine that the encounters between hydroxylases and oxygen occur at cell boundary. In such a way, hydroxylases can be ideally pictured as antennae, localized on the membrane wall, chasing for the (outer) oxygen molecules.

From hereon we will simply refer to VEGF-A as to VEGF.

Importantly, the oxygen is initially imagined to fill the external reservoir, and solely reintegrated via the VEGF pathway. This is clearly a simplifying assumption, as in general one should also account for the oxygen supply via blood flux. On the other hand, we are here interested in the ability of the cell to self-organize when exposed to stressing condition, as it is the lack of oxygen income, a working hypothesis which justifies the assumed scenario.

In deriving eqs. (18) correlations have been neglected, a legitimate choice in the mean-field approximation and whose validity has been tested versus numerical simulation in a number of models that follow the same conceptual scheme[START_REF] Mckane | Amplified biochemical oscillations in cellular systems[END_REF][START_REF] Di Patti | Can a microscopic stochastic model explain the emergence of pain cycles inpatients?[END_REF].

In the transient dynamics, HIF-α molecules are indeed present and translate into VEGF entities, upon combination with HIF-β. VEGF molecules can remain temporarily active also when the HIF-α has been eliminated, and so drive a consequent oxygen influx. There is hence a certain degree of inertia which increments the initial oxygen amount. The effect becomes more pronounced as gα gets larger, more HIF-α populating the cell during the transient dynamics, before degradation has occurred. However, provided gα is small enough, the cell reaches a steady state configuration, with no residual amount of HIF-α molecules. For this reason, we refer to this configuration as to the normoxia regime.


Frank J Poelwijk

Sorin Tȃnase-Nicola

Daniel J Kiviet

Sander J Tans
email: tans@amolf.nl.

Sorin Tǎnase-Nicola

Reciprocal sign epistasis is a necessary condition for multi-peaked fitness landscapes

Having multiple peaks within fitness landscapes critically affects the course of evolution, but whether their presence imposes specific requirements at the level of genetic interactions remains unestablished. Here we show that to exhibit multiple fitness peaks, a biological system must contain reciprocal sign epistatic interactions, which are defined as genetic changes that are separately unfavorable but jointly advantageous. Using Morse theory, we argue that it is impossible to formulate a sufficient condition for multiple peaks in terms of local genetic interactions. These finding indicate that systems incapable of reciprocal sign epistasis will always possess a single fitness peak. However, reciprocal sign epistasis should be pervasive in nature as it is a logical consequence of specificity in molecular interactions. The results thus predict that specific molecular interactions may yield multiple fitness peaks, which can be tested experimentally.

Introduction.

Organisms are highly integrated functional systems, where change in one component often affects the functioning of other components and the system as a whole.

That this interdependence of parts presents challenges for evolutionary processes was already remarked by Darwin [START_REF] Darwin | On the Origin of Species by Means of Natural Selection[END_REF]. Manifestations of these dependencies at the genetic and physiological level are pleiotropy, where a single gene or mutation has an effect on multiple phenotypic traits, and epistasis, which classifies the interaction between genes or mutations in their effect on fitness. The molecular basis for epistatic effects often lies in the physical interactions within or between gene products [START_REF] Kogenaru | Revealing evolutionary pathways by fitness landscape reconstruction[END_REF]), but it can for example also result from the interplay between protein stability and catalytic activity [START_REF] Depristo | Missense meanderings in sequence space: a biophysical view of protein evolution[END_REF], from Watson-Crick base pairing within an RNA molecule [START_REF] Kirby | Maintenance of pre-mRNA secondary structure by epistatic selection[END_REF], or from interactions between modular metabolic networks [START_REF] Segre | Modular epistasis in yeast metabolism[END_REF].

Epistatic interactions between parts may constrain the evolutionary optimization of a biological system. For instance, they may limit the number of trajectories towards a fitness optimum that are selectively accessible, with the latter meaning that each mutational step along the trajectory must be fitness increasing [START_REF] Weinreich | Perspective: Sign epistasis and genetic constraint on evolutionary trajectories[END_REF]. It has been shown that a necessary and sufficient condition for such a reduced accessibility is the presence of so-called sign-epistatic interactions [START_REF] Weinreich | Perspective: Sign epistasis and genetic constraint on evolutionary trajectories[END_REF]. Signepistasis refers to the situation where the fitness effect of a mutation A→a can have a different sign depending on the genetic background B or b (Fig. 1, left). As a consequence, the mutations must occur in a specific order for the paths to be selectively accessible. Sign-epistatic interactions between mutations have been observed empirically, by reconstructing neighboring genotypes and determining their fitness [START_REF] Lunzer | The biochemical architecture of an ancient adaptive landscape[END_REF][START_REF] Bridgham | Evolution of hormonereceptor complexity by molecular exploitation[END_REF][START_REF] Weinreich | Darwinian evolution can follow only very few mutational paths to fitter proteins[END_REF].

A different situation arises when the fitness landscape contains multiple peaks. By definition, not a single path is selectively accessible by single-mutation steps between two peaks (i.e., all paths exhibit fitness decreases). Escape from a suboptimal peak may still be possible, but it involves significantly increased waiting times as specific combinations of genetic changes must appear in a population [START_REF] Phillips | Waiting for a compensatory mutation: phase zero of the shiftingbalance process[END_REF][START_REF] Weissman | The rate at which asexual populations cross fitness valleys[END_REF]), or requires large population sizes that support many mutant genotypes at a given time [START_REF] Weinreich | Perspective: Sign epistasis and genetic constraint on evolutionary trajectories[END_REF]. Entrapment on sub-optima and its relation to landscape ruggedness has principally been studied theoretically, for instance within the mathematical framework of NK landscape models that were pioneered by Kauffman [START_REF] Kauffman | The Origins of Order: Self-organization and Selection in Evolution[END_REF]. Empirical evidence for multiple adaptive peaks is scarce, though indications have been found in experimental evolution studies [START_REF] Korona | Evidence for multiple adaptive peaks from populations of bacteria evolving in a structured habitat[END_REF][START_REF] Salverda | On the natural and laboratory evolution of an antibiotic gene[END_REF]), and by genotype-phenotype mapping [START_REF] Poelwijk | Evolutionary potential of a duplicated repressor-operator pair: simulating pathways using mutation data[END_REF][START_REF] Poelwijk | Empirical fitness landscapes reveal accessible evolutionary paths[END_REF][START_REF] De Visser | Exploring the Effect of Sex on Empirical Fitness Landscapes[END_REF][START_REF] Lozovsky | Stepwise acquisition of pyrimethamine resistance in the malaria parasite[END_REF].

In this paper we investigate the relation between epistasis and multiple peaks. It has long been recognized that interaction between mutations can cause landscape ruggedness [START_REF] Wright | The roles of mutation, inbreeding, crossbreeding and selection in evolution[END_REF], and the involvement of epistasis in multiple peaks has been discussed [START_REF] Weinreich | Perspective: Sign epistasis and genetic constraint on evolutionary trajectories[END_REF]. Here we show that when the landscape harbors multiple peaks, it necessarily contains an epistatic motif that is referred to as reciprocal sign epistasis (Fig. 1, right) [START_REF] Poelwijk | Empirical fitness landscapes reveal accessible evolutionary paths[END_REF]. Put differently, multiple fitness peaks cannot be formed by regular sign-epistasis or magnitude epistasis alone. Reciprocal sign epistasis occurs when two mutations are separately deleterious but jointly enhance fitness: mutation A→a can be advantageous or deleterious, depending on the genetic background B or b, while conversely mutation B→b can also be advantageous or deleterious, depending on the genetic background A or a. We refer to the appendix for a more formal definition of reciprocal sign epistasis. An intuitive physical model for a reciprocal sign-epistatic interaction is that of a lock and a key: if one changes either the lock or key, the key will no longer fit, but changing them both can produce a new matching lock and key combination. Reciprocal sign epistasis thus implies the existence of underlying functional and structural elements within the biological system that are distinct from those that give rise to regular sign-epistasis. Specific lock-key recognition capabilities are pervasive in organisms, and we therefore anticipate that reciprocal signepistasis is a crucial property of biological systems.

Results

To demonstrate that reciprocal sign epistasis is a necessary condition for the existence of multiple peaks, we have to show that multiple peaks imply the existence of a reciprocal sign epistasis motif somewhere in the landscape. Here we describe a graphical method to determine these motifs, of which there may be multiple within the landscape, by considering direct paths from one maximum P to another maximum Q, in a bi-allelic multi-locus system. We consider a path between P and Q to be direct when the number of single mutations along it is equal to the Hamming distance d H , the number of sequence positions where P and Q differ. In the appendix we provide a more rigorous version of this proof.

If we take a direct path between two maxima in a bi-allelic multi-locus system (Fig. 2, solid line), we will find the fitness minimum (indicated with I in Fig. 2) at some position, defined by the number of mutations away from the first peak. We denote the mutation leading into this minimum by A→a, and the following mutation leading away from it by B→b. We can now try to optimize the minimum fitness of this path by exchanging these two mutations, so B→b first, followed by A→a. In many cases the reversal will cause the minimum to shift to a different position along the path between P and Q (e.g. to position II in Fig. 2). If it does, another mutation (C→c) will now lead into or out of the new minimum, and one can thus perform another reversal. This optimization process can be continued until the minimum remains located at the same position (position III in Fig. 2). At this point, no new reversal can be applied around the minimum, and the optimization procedure thus comes to an end. As the graph illustrates, these two mutations exhibit reciprocal sign epistasis by definition, thus proving the presence of this motif.

In the above method, different starting points may locate different reciprocal sign epistasis motifs that can exist within the landscape. To establish the presence of one specific reciprocal sign epistasis motif, namely the one that is highest in fitness, one can follow a simpler argument. Among the paths connecting two peaks, there always exists one unique path that has the highest fitness minimum. When the two mutations around this minimum are exchanged, the new minimum must be lower. Any other outcome would contradict with the assumption that the unique starting path has the highest fitness minimum. Hence reciprocal sign epistasis is established.

We also explored whether the presence of reciprocal sign epistasis implies the existence of multiple peaks, and more generally, whether a sufficient criterion for multiple peaks can be formulated in terms of local two-way genetic interactions. For this purpose we first consider the simple case of a two-allele system detailed in Fig. 3, where each axis represents one locus, the two points along each axis represent the two alleles, and the size of the points denotes the fitness. The two loci indicated by the grayed front face exhibit reciprocal sign epistasis, with the high fitness bottom-left and double mutant top-right separated by low-fitness single mutants. However, a mutation in the third locus does provide an accessible path to higher fitness from the bottom-left, indicating that it does not represent a peak in this landscape, which thus has just a single peak. This simple example shows that reciprocal sign epistasis does not imply the presence of multiple peaks.

The question remains whether other criteria can be formulated that indicate the presence of multiple peaks in terms of local two-way interactions. The generic possibility of having bridging paths along unexplored loci suggests that this may be difficult, and that only a global mapping of the relevant genotype space with all its higher-order genetic interactions provides a sufficient criterion. Insight into the possibility of defining a sufficient criterion can be found using Morse theory [START_REF] Milnor | Morse Theory[END_REF]. Morse theory is a powerful tool to investigate the topology of a smooth manifold (surface or space) using the properties of differentiable functions on the manifold. In particular it provides a set of inequalities relating the number of critical points of the function (minima, maxima, saddles etc). For example, for functions defined on the continuous real manifold of N dimensions, a maximum is a critical point with N negative directions and zero positive directions. Morse theory indicates that the number of saddles with only one positive direction is larger than the number of maxima minus one. This inequality can be compared to our proof (Fig. 2 and appendix) that at least one reciprocal sign epistasis motif, which has similar properties as a saddle, must exist when there are two maxima.

More generally, Morse inequalities state that the number of maxima is larger than the number of saddles of order one (saddle with one positive direction) minus the number of saddles of order two (saddles with two positive directions). For instance, if there would be N saddles of order 1, then one could still have just one peak, as there could be also N saddles of order 2. This statement has a straightforward consequence, namely that to know the number of maxima one has to know all the saddles of order two, and hence have information on the complete surface. By analogy to fitness landscapes, this statement would imply that to determine the number of fitness peaks, one must have exhaustive information on the complete landscape, again indicating that information on local two-way interactions is not sufficient. This argumentation is intuitively captured in the illustration in Fig. 3: even if reciprocal sign epistasis exists, other parts of the landscape can bridge high points and hence smooth-out the ruggedness.

We are not aware of a generalization of Morse theory that applies directly to genotype spaces. Morse theory has been extended to discrete functions defined on spaces obtained by gluing together points, lines, triangles and larger dimensional counterparts [START_REF] Forman | Morse theory for cell complexes[END_REF], which represents a step towards this aim. It would be of great interest to formally extend Morse theory to genotype spaces. Taken together, the arguments suggest that it is not possible to derive a lower bound on the number of maxima based only on the presence or absence of a small number of local two-way interactions, without knowing higher order genetic interactions.

Discussion

The presence of multiple peaks may be considered to be one of the features within fitness landscapes that most critically affect the course of evolution since it precludes adaptation by stepwise mutational optimization. In multi-peaked landscapes, evolutionary trajectories can become trapped on sub-optimal peaks, which may lead to (temporary) evolutionary stasis. While empirical information on the ruggedness of fitness landscapes remains scarce, significant progress has recently been made by systematically reconstructing neighboring genotypes that affect a specific molecular function, such as the resistance against an antibiotic [START_REF] Burch | Evolution by small steps and rugged landscapes in the RNA virus phi6[END_REF][START_REF] Buckling | Adaptation limits diversification of experimental bacterial populations[END_REF][START_REF] Lunzer | The biochemical architecture of an ancient adaptive landscape[END_REF][START_REF] Bridgham | Evolution of hormonereceptor complexity by molecular exploitation[END_REF][START_REF] Miller | Direct demonstration of an adaptive constraint[END_REF][START_REF] Weinreich | Darwinian evolution can follow only very few mutational paths to fitter proteins[END_REF][START_REF] Fernandez | Fitness landscape of human immunodeficiency virus type 1 protease quasispecies[END_REF]), which has yielded some indications for the presence of multiple peaks [START_REF] Korona | Evidence for multiple adaptive peaks from populations of bacteria evolving in a structured habitat[END_REF][START_REF] Poelwijk | Evolutionary potential of a duplicated repressor-operator pair: simulating pathways using mutation data[END_REF][START_REF] Poelwijk | Empirical fitness landscapes reveal accessible evolutionary paths[END_REF][START_REF] Salverda | On the natural and laboratory evolution of an antibiotic gene[END_REF]. However, conclusively demonstrating that two adaptive peaks are distinct -and not interconnected by a ridge that circumvents the adaptive valley-remains difficult if not impossible, due to the high-dimensional nature of sequence space [START_REF] Whitlock | Multiple Fitness Peaks and Epistasis[END_REF].

Moreover, global fitness maps do not directly provide insight into the elementary properties of the biological system that give rise to specific landscape features. Indeed, the extent to which system architecture and shape of the landscape are fundamentally related remains incompletely answered.

Our current investigation of the relation between the presence of multiple fitness peaks and epistatic interactions has been motivated by these unresolved issues. We find that the presence of multiple fitness peaks requires a specific form of epistatic interactions, where two mutations are separately unfavorable but jointly advantageous. This so-called reciprocal sign epistasis motif [START_REF] Korona | Evidence for multiple adaptive peaks from populations of bacteria evolving in a structured habitat[END_REF][START_REF] Poelwijk | Evolutionary potential of a duplicated repressor-operator pair: simulating pathways using mutation data[END_REF][START_REF] Poelwijk | Empirical fitness landscapes reveal accessible evolutionary paths[END_REF][START_REF] Salverda | On the natural and laboratory evolution of an antibiotic gene[END_REF] reflects a mutual dependence between parts at the functional or physico-chemical level that resembles the specific interaction between a lock and a key. Thus, our results imply that to escape from a suboptimal adaptive peak, selection must overcome at least one key-lock type interaction.

Our analysis also shows that the reverse is not true: the occurrence of reciprocal sign epistasis does not guarantee multiple fitness peaks. Moreover, using Morse theory we show that one cannot formulate a local condition in terms of epistatic interactions for the existence of multiple fitness peaks. In other words, local landscape topology (epistatic interactions) cannot conclusively inform about global landscape properties (the presence of multiple peaks). This means that rigorously proving the existence of multiple peaks requires a systematic mapping of genotype space. While this approach may be unfeasible for most biological functions as sequence spaces rapidly become too large for exhaustive experimentation, it may be achievable for specific biological functions that are confined to a small physical region, and hence entail a limited sequence space. Our findings suggest that multi-peaked fitness landscapes may be observed for biological functions involving specific interactions, such as the binding between transcription factor and DNA binding site or residue interactions within a protein fold, which can be tested experimentally.

Appendix:

In this appendix we prove the following theorem:

In a N-allelic L locus system, reciprocal sign epistasis is a necessary condition for the existence of multiple peaks in the fitness landscape.

Let us first define the different types of epistasis referred to in the main text.

Definition 1: Epistasis: means that the fitness effect of a mutation is conditional on the presence of other mutations (the 'genetic background').

AB aB Ab ab w w → → Δ ≠ Δ
, where w Δ is the fitness difference between two mutational states, and B and b can be considered to be the genetic background for states A and a .

For one mutation in different backgrounds, several classes of epistasis can be discerned:

Definition 2: Magnitude epistasis: means that the magnitude of the fitness effect of a mutation is conditional on the presence of other mutations. 3. A bi-allelic L locus system may contain more than two peaks. However, we may consider the sub-space between any two peaks in the full landscape, defined by the loci where their sequences differ. If this sub-space should contains reciprocal sign epistasis, then the full landscape also contains it. Hence, we can restrict ourself to two peaks.

So finally, we are left to prove the following theorem:

In a bi-allelic L locus system, reciprocal sign epistasis is a necessary condition for the existence of two peaks at a distance L in the fitness landscape.

The following proof will show that along a path containing the highest fitness minimum in the landscape, at least two mutations exhibit reciprocal sign epistasis. Because all paths between two maxima necessarily go through a minimum, the landscape must contain reciprocal sign epistasis.

1) There are two maxima at distance L in a biallelic L locus system.

2) We consider a random direct path (length L) from one maximum to the other one.

3) This path necessarily contains a minimum.

  Sign epistasis: means that the sign of the fitness effect of a mutation is conditional on the presence of other mutations. Reciprocal sign epistasis: means that the sign of the fitness effect of mutation A→a is conditional on whether the state of another locus is b or B , We start by narrowing down the scope of systems for which we need to prove the theorem:1. If reciprocal sign epistasis must occur for direct paths between two peaks, then reciprocal sign epistasis must also occur for all paths, as they include both direct and indirect paths. Hence, we can restrict ourself to direct paths between two peaks.2.The set of direct paths with length d between two points in a N-allelic L locus system is equal to the set of direct paths of length d in the bi-allelic L locus system. This is rather trivial: direct paths between between two points in sequence space ( 0 = H d d -) require only one substitution at each locus, thus only two alleles. Hence we can restrict ourself to the bi-allelic case.

4)FiguresFigure 1 :

 1 Figures

Figure 2 :

 2 Figure 2: Procedure for finding the reciprocal sign epistasis motif. A random direct path

Figure 3 :

 3 Figure 3: Schematic diagram illustrating that reciprocal sign epistasis is not a sufficient

Acknowledgements

This work is part of the research program of the Stichting voor Fundamenteel Onderzoek der Materie, which is financially supported by the Nederlandse Organisatie voor Wetenschappelijke Onderzoek.


Anthony D Fox

Roy King

Moulting Pochard

Aythya Ferina

Tufted Duck

A Fuligula At

Body mass loss amongst

Keywords: wing moult, moult migration, phenotypic plasticity, fat stores, energy budget Zusammenfassung Verlust von Körpermasse bei mausernden Tafel-(Aythya ferina) und Reiherenten Anforderungen während der mauserbedingten Flugunfähigkeit bei Tafel-und

We assessed the contribution of endogenous fat stores to meeting energetic needs during the flightless moult period in Pochard and Tufted Duck, by regressing mass on the progression of wing moult, measured by primary length. Pochard lost between 22.1% (males) and 24.2% (females) of body mass and female Tufted Duck 12.2% during wing moult at Abberton Reservoir, Essex. Based on a 27 day flightless period, Pochard lost on average 8.3-8.4 g per day and Tufted Duck 3.2-3.4 g per day, assumed to be the result of fat expenditure. Assuming daily energy expenditure (DEE) to be 1.7x basal metabolic rate (BMR), the contribution from burning endogenous fat equated to 37-40% of Pochard and 19-20% of Tufted Duck DEE during remiges moult. One frequently re-trapped female Pochard expended fat stores that almost fulfilled her entire energetic requirements during moult if she rested and did not feed for most of the moult. These results confirm that Pochard and Tufted Duck accumulated fat stores that helped meet energetic needs during the flightless wing moult period. This endogenous source of energy may free them to exploit habitats that are safe from predation but which may not fully supply the energetic needs of moulting ducks from exogenous sources.

Reiherenten durch Regression der Masse gegen das Fortschreiten der Flügelmauser, gemessen an der Länge der Handschwingen. Am Abberton Reservoir in Essex verloren Tafelenten zwischen 22,1% (Männchen) and 24,2% (Weibchen) Körpermasse, während Reiherentenweibchen 12,2% Körpermasse verloren. Gemessen an einer flugunfähigen Periode von 27 Tagen verloren Tafelenten durchschnittlich 8,3 -8,4 g pro Tag und Reiherenten 3,2 -3,4 g pro Tag, was vermutlich auf die Nutzung von Fettdepots zurückzuführen ist. Unter der Annahme eines täglichen Energiebedarfs (DEE) von 1,7x Grundstoffwechselrate (BMR) beträgt der Verbrauch von endogenen Fettreserven während der Mauser der Flugfedern 37-40 % des DEE bei Tafelenten und 19-20 % bei Reiherenten. Eine mehrfach wieder gefangene weibliche Tafelente verbrauchte Fettreserven, die allein nahezu ihren gesamten Energiebedarf während der Mauser decken könnten, sofern sie rastet und nicht nach Futter sucht. Diese Ergebnisse bestätigen, dass Tafel-und Reiherenten Fettreserven anlegen, die ihnen helfen, die flugunfähige Zeit der Mauser zu überbrücken. Diese endogenen Energiereserven erlauben es ihnen Habitate aufzusuchen, die sicher vor Fressfeinden sind, aber nicht ausreichend Futterquellen bieten, um den Energiebedarf während der Mauser zu decken.

Introduction

Most ducks replace wing feathers simultaneously, rendering them temporarily flightless; males and nonbreeding females (sub-adult, other non-or failed breeders) often do so remote from breeding areas (the "moult migration" of [START_REF] Salomonsen | The moult migration[END_REF]). Protein synthesis is energetically costly [START_REF] Portugal | Annual changes in body mass and resting metabolism in captive barnacle geese (Branta leucopsis): the importance of wing moult[END_REF]) and feather protein, keratin, is characterized by its high sulphur-containing amino-acid content, making wing moult a period of potential nutrient and energetic stress [START_REF] Hanson | The dynamics of condition factors in Canada geese and their relation to nutritional stress[END_REF][START_REF] Hohman | Ecology and management of postbreeding waterfowl[END_REF]. Recent studies suggest that waterfowl show remarkable phenotypic plasticity to meet the needs of the flightless wingmoult period, notably resorbing and reconstructing organs and muscle blocks in anticipation of need (rather than a result of training, [START_REF] Fox | Changes in body mass and organ size during wing moult in non-breeding greylag geese Anser anser[END_REF][START_REF] Portugal | Testing the use/disuse hypothesis: pectoral and leg muscle changes in captive barnacle geese Branta leucopsis during wing moult[END_REF], synthesising flight feathers from protein stores accumulated in the body prior to moult [START_REF] Fox | Isotopic evidence for endogenous protein contributions to greylag goose Anser anser flight feathers[END_REF]) and using stored fat reserves, also accumulated pre-moult, to meet energetic needs during flight feather re-growth [START_REF] Fox | Changes in body mass and organ size during wing moult in non-breeding greylag geese Anser anser[END_REF][START_REF] Kahlert | Effects of feeding patterns on body mass loss in moulting greylag geese Anser anser[END_REF]). Most adaptations are shown by waterbirds exposed to high predation risk, so plasticity responses to flightlessness can be interpreted as mechanisms to enable rapid re-growth of feathers in an environment which enables the birds to avoid predation, without complete reliance on exogenous sources of energy and/or nutrients to balance their daily needs. In other words, where predation risk is high (as, for instance, on land for normally terrestrial grazing geese), waterbirds can potentially reduce mortality rates whilst flightless by restricting feeding time, maximising food retention time and using body stores to help meet energy and nutrient needs during flightless moult [START_REF] Fox | Adjustments to nitrogen metabolism during wing moult in greylag geese[END_REF].

It might therefore be expected that waterbirds exploiting very extensive aquatic habitats, characterised by abundant food resources and low predation risk, would exhibit reduced reliance upon endogenous stores, particularly to meet the energetic costs of moult. In moulting seaducks, a combination of an abundant and unrestricted benthic food supply and lack of predators in the open sea would predict no loss of body mass and particularly of fat stores during moult; as is the case for moulting Common Scoter Melanitta nigra in Danish waters [START_REF] Fox | Changes in body mass and organ size during remigial moult in common scoter Melanitta nigra[END_REF]. Similar patterns seems likely in freshwater diving ducks that moult remiges on very large water bodies, especially those that undertake long moult migrations from the breeding areas to moult sites at or near their ultimate winter quarters (e.g. the Aythya ducks of eastern Europe that come to western Europe to moult and winter, [START_REF] Salomonsen | The moult migration[END_REF][START_REF] Fox | Changes in non-breeding distribution and habitat of Pochard Aythya ferina in Britain[END_REF]Wernham et al. 2002). Thus, Tufted Duck A. fuligula and Pochard A. ferina might arrive at moult sites after long migration that would deplete fat stores, making it unlikely that they could invest endogenous energy stores in meeting the needs of feather replacement. However, evidence suggests that at least some Aythya duck species do lose body mass during the wing moult, with 31% depletion of lipid stores in Redheads A. americana (Bailey 1981) and 45% in female Ring-necked Ducks A. collaris [START_REF] Hohman | Annual body weight change in ring-necked ducks (Aythya collaris)[END_REF]), whilst others do not (male Ring-necked Ducks [START_REF] Hohman | Annual body weight change in ring-necked ducks (Aythya collaris)[END_REF], female Lesser Scaup A. affinis Hohman et al. 1992and Canvasbacks A. valisineria Thompson and Drobney 1996[START_REF] Thompson | Diet and nutrition of male canvasbacks during post-reproductive molts[END_REF].

So, do moult migrant Tufted Duck and Pochard in western Europe lose mass due to depletion of fat stores during moult to assist with meeting energy requirements? In situations with abundant food supply and limited predation risk (such as large lakes which form the winter quarters of far higher densities of the same species later in the season), there should be no need to accumulate fat stores for exploitation during the flightless period. Here we test this hypothesis by reporting changes in body mass amongst these two species from birds caught during flightless moult at the large (>500 ha) Abberton Reservoir (51° 49' N, 00° 52' E) an important wintering site for both species in Essex, south east England, which lacks aerial or aquatic predators of ducks.

METHODS

Measurements from live-trapped free-living ducks

Tufted Duck and Pochard were caught opportunistically in baited traps at Abberton during the late summer of 1979. Data are presented here from 27 female (21 flightless, 6 completed moult and able to fly) and 18 males (13 flightless 5 flying) Pochard and 46 female (36 flightless and 10 flying) and 8 flightless male Tufted Ducks. Birds were aged and sexed, fitted with BTO metal rings, the flattened and straighten wing length measured to the nearest mm and body mass determined to the nearest 5 g using a balance by RK.

Notes were taken of the state of moult, especially to differentiate old full grown from new feathers. We reasoned that although the rate of wing feather regrowth may not be linear with time, the length of the longest growing primary gave the best measure of the moult stage of an individual. One female Pochard GK65638, was caught and recaptured 18 times, enabling assessment of feather growth rate and mass change in a single individual throughout the flightless period. We assessed the rate of mass change during wing moult in these two species using linear and quadratic regression.

Modelling mass loss

Previous studies of flightless Greylag Geese Anser anser showed that mass in late moult was almost entirely explained by loos of fat accumulated prior to moult [START_REF] Fox | Changes in body mass and organ size during wing moult in non-breeding greylag geese Anser anser[END_REF]. For the purposes of this analysis, we assumed that mass loss during moult related to fat used to help meet the energetic demands during moult [START_REF] Kahlert | Effects of feeding patterns on body mass loss in moulting greylag geese Anser anser[END_REF]. We have assessed the relative contribution of fat to the overall energy expenditure of the two species during their flightless periods (here considered to be 27 days based on observations). The daily energy contribution of fat during the flightless period was calculated by dividing the mass lost (known precisely in the case of GK65638, modelled by regression for other birds) by 27 days and converted to daily energy by multiplying by 38.9, the energy density of fat in kJ g -1 [START_REF] Pennycuick | Modelling the flying bird[END_REF]. Daily energy expenditure (DEE measured in kJ d -1 ) in Anatidae has been widely estimated as 2.6x basal metabolic rate (BMR) during other parts of the annual cycle [START_REF] Drent | Balancing the energy budgets of arctic-breeding geese throughout the annual cycle: a progress report[END_REF].

This conversion factor incorporates periods of flight, an activity which is energetically costly relative to other activities (10-12x BMR, e.g. [START_REF] Mooij | Behaviour and energy budget of wintering geese in the Lower Rhine area of North Rhine-Westphalia, Germany[END_REF]) but which self-evidently forms no part of the flightless moult period of the ducks. Thus, DEE should be lower than 2.6x BMR this during moult. [START_REF] Bevan | Daily energy expenditure of Tufted Ducks: a comparison between indirect colorimetry, doubly labelled water and heart rate[END_REF] and [START_REF] Leeuw | Wintering Tufted Ducks Aythya fuligula diving for zebra mussels Dreissena polymorpha balance feeding costs within narrow margins of their energy budget[END_REF] experimentally found that DEE of Tufted Ducks amounted to between 1.7 and 2.5x BMR in summer (based on tank experiments indoors and outdoors respectively), which included diving activity, which is also energetically costly. [START_REF] Kahlert | Effects of feeding patterns on body mass loss in moulting greylag geese Anser anser[END_REF] suggested 1.7x BMR was more realistic based on observations of free-living flightless moulting Greylag Geese travelling up to 6 km per day to feed. [START_REF] Portugal | Do captive waterfowl alter their behaviour patterns during their flightless period of moult?[END_REF] have shown that Common Eiders Somateria mollissima (a maritime diving duck) and Garganey Anas querquedula show significant reductions in feeding and locomotion and increases in resting during moult, therefore it seems likely that these two Aythya ducks may also further reduce energy expenditure by resting more and diving less during the flightless period, although this will be offset by the additional energetic costs of feather synthesis [START_REF] Portugal | Annual changes in body mass and resting metabolism in captive barnacle geese (Branta leucopsis): the importance of wing moult[END_REF]). We compare the calculated DEE of the two species during moult based on 2.6x, 1.7x and 1.0x BMR to assess the degree to which Tufted Duck and Pochard could meet their DEE under the differing scenarios, varying from normal energy consumption and to highly elevated levels of resting, to no feeding (i.e. no diving).

RESULTS

Body mass (M) significantly declined with wing length (w) amongst female and male Pochard caught while flightless (Fig. 1). The best fit regression was M = 1164.8 -2.53w (F 1,20 = 12.4, r = 0.63, P = 0.002) for females and M = 1236.70 -2.55w (F 1,11 = 7.3, r = 0.65, P = 0.02) for males. Taking the moult period as being that between the shortest and longest wing measurements on flightless birds and substituting these values in the regression models equates to losses of 227 g and 224 g in females and males respectively, or 24.2% and 22.1% of the original weights at the start of moult. In the case of GK56538, the regression model had the formula M = 1198.9 -3.06w (F 1,18 = 7.3, r = 0.97, P << 0.0001).

There was a significant decline in body mass with increasing wing length in female Tufted Ducks caught during flightless moult (M = 798.1 -0.97w F 1,36 = 6.6, r = 0.40, P = 0.014) and a similar rate of decline among males, although the male regression was not significant due to low sample size (M = 815.1 -0.93w (F 1,7 = 2.2, r = 0.52, P = 0.184, Fig. 1). Using the same approach as for Pochard, this equates to losses of 86.9 g in females and 92.9 g in males, 12.2% and 12.6% of the original weights at the start of moult.

If all mass loss in moulting Pochard reflects fat use, the amount of energy available to both sexes from fat stores during moult amounted to c.25% of their total daily requirements at normal levels of DEE, and rises to c.38% if energy costs were 1.7x those of BMR (Table 1). If female GK65638 had ceased foraging and simply subsisted at normal BMR for the duration of the flightless period, her fat stores could have met 92% of her energy needs during her flightless period (Table 1).

Based on the same assumptions as for Pochard, moulting Tufted Duck of both sexes were able to meet only c.12% of their total daily requirements based on 2.6x BMR, raising to c.19% if energy costs were 1.7x those of BMR and just under one third of DEE if energetic costs were those of resting metabolism (Table 1). synthesis (1.7x BMR). In this species, endogenous fat stores could supplement, but not completely meet, energy demands during flightless moult. In contrast, the "trap happy" Pochard GK65638 was potentially able to meet over 90% of her energy requirements at resting metabolism from fat storage alone. We should be extremely cautious in the interpretation of patterns shown by this one individual, because her frequent recapture could suggest poor body condition or otherwise aberrant behaviour. Although she benefited from bait at the trap, her case suggests that ducks could meet their DEE for flight feather replacement without feeding. Our data likely under represent the extent of body mass loss during the flightless period, since few birds (not even GK65638) were caught immediately after shedding flight feathers, when body mass is likely to be at its highest point during the moult. Furthermore, the flightless period may be shorter than the 27 days we used to calculate DEE (3-4 weeks reported in Cramp and Simmons 1977, but Mallard Anas platyrhynchos can vary between 22 and 29 days, [START_REF] Panek | Remex growth and body mass of Mallards during wing molt[END_REF]. If the flightless period was as short as 22 days, GK65638 could have met all her energy needs at DEE of 1x BMR with substantial unused fat stores, and the "average" Pochard could have met c.80% of their DEE from endogenous fat stores.

Species and sex differences in mass loss

There was no evidence of differences in mass loss between the sexes within the two species considered here at the same moulting site, although the small sample size amongst Tufted Duck restricts our ability to make firm conclusions for this species. Pochard consumed a quarter and Tufted Duck one eighth of their pre-moult body mass, suggesting consistent between-species differences at Abberton Reservoir. Why the difference? Typically, both of these bird species take advantage of feeding by night on areas that are often subject to disturbance by day, for example, waters subject to fishing [START_REF] Bell | The game-fishing season and its effects on overwintering wildfowl[END_REF], intensive water sports [START_REF] Fox | Food supply and the effects of recreational disturbance on the distribution and abundance of wintering Pochard on a gravel pit complex in southern Britain[END_REF], shooting disturbance [START_REF] Evans | Does shooting disturbance affect diving ducks wintering on large lakes? A case study on Lough Neagh, Northern Ireland[END_REF] or building and redevelopment [START_REF] Marsden | Impact of Disturbance on Waterfowl Wintering in a UK Dockland Redevelopment Area[END_REF]. Hence, in western Europe these species may exploit risky, but food-rich habitats outside of the breeding season that are generally subject to human disturbance or sources of high predation risk, and use safer, geographically separated habitats for roosting and resting at other times during the 24 hour cycle [START_REF] Fox | Food supply and the effects of recreational disturbance on the distribution and abundance of wintering Pochard on a gravel pit complex in southern Britain[END_REF]. This may mean that, in order to exploit such habitats during the vulnerable flightless moult period, evolution has favoured a strategy in these species that energetically buffers birds against the costs of flight feather replacement and enables them to exploit habitats which may not fulfil all of their energy requirements but which do offer safety from predation. But why the difference between these two species and other Aythya species, such as the Canvasback? Perhaps the explanation is diet, because Tufted Duck consume more animal prey with higher fat and protein content than the more omnivorous Pochard. Pochard tend to occur in larger groups than Tufted Ducks during the non-breeding period and that these aggregations more often feed and roost on different water bodies than Tufted Duck. If Tufted Ducks are subject to less predation when flightless than are Pochard, then the interspecific differences may reflect the ability of their moulting sites to provide adequate energy and protein to sustain moulting ducks without risk of predation. Whatever the explanation, these Anatid species show remarkable plasticity in their ability to renew flight feathers simultaneously without being forced to use feeding habitats that may expose them to high predation risk, but we still need to better understand the ecophysiology of Anatids during the critical wing moult period. 1. Estimates of daily energy consumption of moulting Pochard at Abberton Reservoir, Essex, southeast England. Body mass changes are shown for a single female GK65638 which was re-trapped almost daily throughout the flightless moult period, as well as for "average" females and males based on the regression models described in the text. It is assumed that the body mass loss equates entirely to fat consumption and that 1 g fat provides 38.9 kJ energy. Predicted Basal Metabolic Rate (measured in kJ per day) is based on the equation BMR = 3.06 x mass 0.74 and daily energy expenditure is based on different multiples of BMR as indicated, expressing the potential daily source of energy from fat stores as a percentage of daily energy needs.

Figure 1 .

 1 FIGURE and CAPTION

TABLE Total

 Total

			Imputed daily	Energy	Initial	Predicted daily energy
		imputed fat	fat mass	supplied	body	consumption based on initial
		mass	consumption	by this	mass	body mass estimated by
		consumption	based on a	amount	(g)	different levels of energy
		(g)	flightless period	of fat (kJ		consumption relative to BMR
			of 27 days	per day)		(percent of total energy needs)
			(g per day)
						2.6 x	1.7 x	1.0 x
						BMR	BMR	BMR
	Pochard	310	11.5	446.7	940	1261.3	824.7	485.1
	Female					(35.4%)	(54.2%)	(92.1%)
	GK65638
	Average	228	8.4	327.6	945	1266.2	827.9	487.0
	Pochard					(25.9%)	(39.6%)	(67.3%)
	Female
	Average	224	8.3	322.2	1014	1333.6	872.0	512.9
	Pochard					(24.2%)	(37.0%)	(62.8%)
	Male
	Average	87	3.2	125.2	711	1026.1	670.9	394.6
	Tufted Duck					(12.2%)	(18.7%)	(31.7%)
	Female
	Average	93	3.4	133.9	737	1053.6	688.9	405.2
	Tufted Duck					(12.7%)	(19.5%)	(33.0%)
	Male

Table

ACKNOWLEDGEMENTS

Thanks to Johnny Kahlert and two anonymous referees for comments and improvements on an earlier version and to the then Essex Water for permission to WWT to catch and ring ducks at Abberton Reservoir.

DISCUSSION

Do Pochard and Tufted Duck lose mass during moult? [START_REF] Thompson | Nutritional consequences of molt in male canvasbacks: variation in nutrient reserves and digestive tract morphology[END_REF] showed that Nearctic male Canvasbacks gradually catabolised somatic lipids from arrival on moulting sites through their mid-remigial moult. They argued that if Canvasbacks relied on lipid stores to meet the energy requirements of moult, they should have stored fat after their arrival at moulting sites and prior to the flightless period; however, there was no evidence for this fattening schedule [START_REF] Thompson | Intensity and chronology of post-reproductive molts in male Canvasbacks[END_REF]. Instead, lipogenesis began during mid-remigial growth, suggesting that male Canvasbacks primarily rely on their diet to meet energetic needs during moult. In contrast, the results presented here suggest more dramatic and consistent mass loss amongst two Palearctic Aythya species throughout the period of flightlessness. Pochard lost about 23% of their initial biomass during the moult period and Tufted Duck lost 12%. We infer that this mass loss reflects use of stored fat for energy, because changes in breast and leg muscles in Canvasbacks could not explain such mass loss [START_REF] Thompson | The nutritional ecology of molting male Canvasbacks (Aythya valisineria) in central Alberta[END_REF]. Evidence from other moulting waterfowl suggests such rapid changes in mass are not linked to changes in protein, water or other body components (e.g. Young and Boag 1982;[START_REF] Thompson | Nutritional consequences of molt in male canvasbacks: variation in nutrient reserves and digestive tract morphology[END_REF] but are related to fat metabolism [START_REF] Fox | Changes in body mass and organ size during wing moult in non-breeding greylag geese Anser anser[END_REF]) despite the caution of [START_REF] Piersma | Pre-migratory "fattening" usually involves more than the deposition of fat alone[END_REF].

The mass of both Pochard and Tufted Duck immediately after completion of remiges moult showed some signs of stabilisation and even increased, suggesting that the mass loss we observed was specifically associated with the flightless period (see Fig. 1). Although we lack large samples of birds caught in the prelude to moult, mass in both species was higher in the pre-moult period than in mid summer or early autumn suggesting a pre-moult acquisition of fat that [START_REF] Thompson | Nutritional consequences of molt in male canvasbacks: variation in nutrient reserves and digestive tract morphology[END_REF] failed to find in the Canvasback. In this respect, Pochard and Tufted Duck show similarities to dabbling ducks [START_REF] Folk | The weight of the Mallard Anas platyrhynchos and its changes in the course of a year[END_REF][START_REF] Hay | Molting biology of male Gadwall at Delta, Manitoba[END_REF]Young and Boag 1982;[START_REF] Pehrsson | Effects of body condition on molting Mallards[END_REF][START_REF] Sjöberg | The flightless period of free-living male Teal Anas crecca in northern Sweden[END_REF][START_REF] Panek | Remex growth and body mass of Mallards during wing molt[END_REF] and southern moulting geese [START_REF] Fox | Changes in body mass and organ size during wing moult in non-breeding greylag geese Anser anser[END_REF][START_REF] Portugal | Annual changes in body mass and resting metabolism in captive barnacle geese (Branta leucopsis): the importance of wing moult[END_REF].

Contribution of fat stores to total energy needs during moult

Table 1 suggests that Tufted Duck could fulfil c.19% of their DEE during moult from fat stores, assuming some energy expenditure for feeding and enhanced energetic consumption to meet the demands of feather