Édouard Challe

IS THE STUDY OF BUSINESS-CYCLE FLUCTUATIONS "SCIENTIFIC?"

Keywords: theory of fluctuations, general equilibrium, fiscal multipliers. 1. See, for example, Reis

The study of macroeconomic fluctuations assumes that the behavior of the whole (aggregates) cannot be reduced to the sum of the parts (agents, markets). This is because interdependencies between markets can substantially amplify, or on the contrary dampen, shocks that at any time disturb the equilibrium. The understanding of general-equilibrium effects, on which direct evidence is limited, which are empirically blurred by multiple potential confounding factors, and for which controlled experiments are almost impossible to design, is necessarily more conjectural than the study of individual behavior or of a specific market. However, ignoring these effects because they do not have the same degree of empirical certainty as a directly observed microeconomic effect can lead to serious policy mistakes.

Busi ness-cycle macroeconomics has been the subject of much criticism in recent years, to the point that it is often perceived from the outside as a field in an irremediable state of crisis 1 . I will focus here on the criticism, or rather the cluster of criticisms, potentially the most destructive, which consists of questioning the mere "scientificity" of business-cycle macroeconomics, not only with regard to other sciences (criticism which, whatever one thinks about it, is old) but in light of the recent evolution of the economics itself, especially as regards its closer relationship to data. There are two main sides to this general criticism:

-The study of macroeconomic fluctuations would not have achieved the "empirical turn" characteristic of mature disciplines, even though examples of such a turning point are numerous in related fields such as labor economics, development economics, as well as in many areas of microeconomics such as industrial organization or corporate finance. Having missed this empirical turn, business-cycle macroeconomics would still consist of speculating on plausible causalities, conjectures, imaginary worlds that are potentially far removed from the one in which we live;

-Moreover, and this is partly a variation of the previous point, the theory of macroeconomic fluctuations would face an almost insurmountable problem of falsification: to the extent that too few data is available to choose among too many macroeconomic models, the stock of available models supposedly accumulates without limit over time without any effective sorting taking place. To borrow Noah Smith's expression, macroeconomists would tend to "cover all the bases" [START_REF] Buchanan | Debating what's wrong with macroeconomics[END_REF], demultiplying models and their associated sets of assumption indefinitely, instead of selecting a small number of relevant models.

These criticisms are severe, but are they truly justified? In any case, they do not seem to take into account an essential dimension of the study of fluctuations, which distinguishes it from other fields of economics: the importance it attaches to the strategic interactions between agents as well as to the general-equilibrium effects that take place across different markets. This is what makes macroeconomics in general, and business-cycle macroeconomics in particular, truly special: it is based on the mere notion that the behavior of the whole (macroeconomic aggregates) cannot be reduced to the sum of the parts (the agents, the markets). This is because the various interdependencies between agents and between markets can substantially amplify, or on the contrary dampen, shocks that at any moment disturb the equilibrium. The understanding of strategic interactions and general equilibrium effects, on which direct evidence is limited, which are empirically blurred by multiple potential confounding factors, and for which controlled experiments are almost impossible to design, is necessarily more conjectural than the study of individual behavior or of a specific market. However, ignoring these interdependences on the grounds that we cannot reach about them the same degree of empirical certainty as about an isolated microeconomic mechanism would not only prevent us from understanding certain complex and large-scale phenomena (such as the "Great Recession"), but can also lead to misguided economic policy recommendations. In what follows I develop these two points by relying on a critical discussion of the recent literature.

Strategic Interactions and General Equilibrium Effects: Between Amplification and Dampening of Aggregate Shocks

In most cases, we do not observe a "macroeconomic shock" that alone can explain the extent of an economic crisis. The Great Depression of the 1930s was preceded by a modest stock market crash, of which no one could have anticipated the effects. The Great Recession that followed the 2008 crisis followed a major financial shock but was quickly contained by the concerted action of the major central banks; this shock alone cannot explain the depth and duration of the Great Recession, even in the United States. More generally, we do not have direct evidence of large shocks along the business cycle that alone could explain its amplitude. If production and employment vary so much over the business cycle, it must be that the economic system contains the seeds of its own instability, by amplifying the impact of small disturbances. Such amplification mechanisms are difficult to identify empirically because they generally involve several mechanisms simultaneously set in motion and generate co-movements of all macroeconomic variables. Understanding such intricacies is usually impossible without a fully specified general-equilibrium model, which explains why business-cycle analysis gives a prevalence to macroeconomic theory over a more inductive, empirical approach. Let me illustrate this point using the main three propagation mechanisms that have be argued to have contributed to the depth and duration of the Great Recession.

The liquidity trap and the deflationary spiral

The liquidity trap is defined as a situation in which the abundance of reserve money within the banking system causes the nominal interest rate on the interbank market to fall to the level of the interest rate on the excess reserve that private banks hold on their account with the central bank. This occurs precisely when the central bank is attempting to implement the maximum level of monetary accommodation, so in a situation where the interest rate on reserves is itself kept close (and potentially slightly below) the rate of return on bank notes, namely zero (abstracting from the cost of storing bank notes). At this point "conventional" monetary policy becomes inoperative and any macroeconomic shock is magnified by the deflationary spiral depicted in Chart 1: a falling aggregate demand depresses output and inflation; at a given nominal interest rate, these deflationary pressures cause a rise in the real interest rate, thereby reinforcing the initial fall in aggregate demand, and so on. This feedback loop has been the subject of a large literature since the pioneering contribution of [START_REF] Krugman | It's Baaack! Japan's slump and the return of the liquidity trap[END_REF] and is today one of the main explanatory model for the depth and duration of the Great Recession in the United States and the euro area. 2 The deflationary spiral depicted in Chart turns out to be particularly difficult to directly measure empirically -far more than the effect of limited disruption in a particular market. This spiral involves a number underlying macroeconomic blocs (the Phillips curve, the Fisher relation, demand-determined output), each of which with its own identification challenges. Given the inherent complexity of the economic mechanism at work, attempts at empirically evaluating this deflationary spiral have essentially adopted one of the following approaches:

-The first approach is to test a specific implication of the propagation mechanism under consideration, which clearly distinguishes High real interest rate

Low output Low current inflation

Low future inflation

Negative shock it from alternative propagation mechanisms. In the present context, the so-called paradox of toil [START_REF] Eggertsson | The paradox of toil[END_REF], according to which negative supply shocks become expansionary in a liquidity trap (because of their inflationary impact), provides the needed crucial experiment. 3 In this spirit, [START_REF] Wieland | Are negative supply shocks expansionary at the zero lower bound?[END_REF] rejects the liquidity trap hypothesis by showing that negative productivity shocks (earthquakes, oil price shocks) are contractionary even at the zero lower bound. In contrast, [START_REF] Datta | Oil, equity, and the zero lower bound[END_REF] find strong co-movements between oil and equity returns at the zero lower bound, which is supportive of the liquidity trap hypothesis. Even if the issue is not yet empirically settled, it remains that the impact of supply shocks at the zero lower bound provides a clean test of the feedback loop described in Chart 1.

-The second approach is to specify a complete general-equilibrium model, in which the deflationary spiral mechanism is present, and then to estimate it empirically (see, for example, [START_REF] Lawrence | Understanding the great recession[END_REF][START_REF] Gust | The empirical implications of the interest-rate lower bound[END_REF]. This approach makes it possible to measure the full causal chain postulated by the theory and potential then to build alternative scenarios ("counterfactuals") which describe how the economy would have behaved if this causal chain had been broken (say, if the bank central could have implement very negative interest rates).

It is clear that in both cases economic theory plays a preponderant role. In the first case, a complete dynamic general equilibrium model is necessary to formulate a testable implication of the considered mechanism; in the second, the full same model (potentially augmented with addition features) is itself estimated on historical data. In any case the deflationary spiral does not spontaneously show up in macroeconomic time-series: it is primarily a theoretical construct and therefore, from the outset, an interpretation of these time-series.

The precautionary-saving feedback loop

A second amplification mechanism, which can play simultaneously or independently of the previous one, involves the precautionary saving behavior of households and the way in which it interacts with unemployment risk over the business cycle. This spiral is summarized in Chart 2. Intuitively, a fall in output that causes employment to fall back raises households' precautionary savings (in anticipation of the increased risk of unemployment); the induced fall in aggregate demand reinforces the initial drop in output and employment, increases the risk of unemployment, and so on. This spiral involves three basic mechanisms. First, output must respond in one way or another to aggregate demand (for example because nominal prices are sticky). Second, labor-market flows and the unemployment risk that they generate must respond endogenously to output changes; this requires a representation of the labor market in terms of worker flows (between employment and unemployment) and not simply in terms of stock (employment). Finally, and perhaps most importantly, households need to be imperfectly insured against the risk of unemployment -otherwise there would be no precautionary motive in the first place and therefore no time-variations in precautionary savings. These three mechanisms are present in various forms, and thus generate the precautionary-saving spiral, in the works of [START_REF] Challe | Precautionary saving and aggregate demand[END_REF], [START_REF] Chamley | When demand creates its own supply: saving traps[END_REF], Den [START_REF] Haan | Unemployment (fears) and deflationary spirals[END_REF], [START_REF] Heathcote | Wealth and volatility[END_REF], [START_REF] Ravn | Job uncertainty and deep recessions[END_REF] and [START_REF] Werning | Incomplete markets and aggregate demand[END_REF], among others.

Empirically measuring the feedback loop depicted in Chart 2 is, again, challenging. Quantitative assessments of the precautionarysaving spiral require from the outset the formulation of a complete dynamic general-equilibrium model in which the three ingredients 2017) calibrate such a model, paying particular attention to the dynamics of the labor market. [START_REF] Challe | Precautionary saving and aggregate demand[END_REF] propose a structural estimation of a related model in order to evaluate the amplifying role of the precautionary motive during the last three recessions in the United States. As far as I know, there is no crucial experiment (of the kind of the paradox of toil in models of the liquidity trap) that would make it possible to directly test the existence of the precautionary-saving spiral: when one departs from the structural estimation of the full general-equilibrium model, only indirect evidence about a particular dimension of the loop is available (relating, for example, to the effect of fluctuations in employment on consumption demand). Thus, just as in the case of the liquidity trap, the precautionary-savings spiral is a plausible propagation mechanism, the amplitude of which can be measured in the data only through the lens of a fully-specified general-equilibrium model. For this reason, the precautionary-saving spiral is best understood a particular way of interpreting the joint dynamics of aggregate demand and unemployment during a recession, which can (and should), be confronted (and possibly associated) with alternative plausible amplification mechanisms.

Note that the work on imperfect insurance and the precautionary motive fully integrates individual heterogeneity into macroeconomic dynamics, recognizing from the outset that different households (in terms of wealth, income, labor market perspective etc.) behave differently, notably in terms of consumption and asset accumulation choices. In particular, a typical result in this literature is that poorer households (that is, households that are close to their own debt constraint) have an individual consumption response to macroeconomic shocks that is stronger than that of richer households, and hence their presence is more likely to set in motion the precautionarysaving feedback loop described in Chart 2. This approach makes it clear that the amplitude of the business cycle and the level of inequalities within a particular economy are fundamentally intertwined. For this reason, this approach allows studying some important economic policy issues that can not otherwise be addressed. Among those issues are the aggregate demand effects of redistributive policies (through taxes or unemployment insurance), which is the focus of the recent work of McKay andReis (2016, 2017).

The credit cycle

A third feedback loop, which is thought to have played an important role in the propagation of the Great Recession, is the so-called "credit cycle" depicted in Chart 3. The modern formulation of the credit cycle dates back to [START_REF] Kiyotaki | Credit Cycles[END_REF]. The theory was later on operationalized into an estimated DSGE model by [START_REF] Iacoviello | House prices, borrowing constraints, and monetary policy in the business cycle[END_REF], who focused on the joint fluctuations of credit and house prices. More recent contributions such as Jeanne and [START_REF] Korinek | Managing credit booms and busts: A Pigouvian taxation approach[END_REF] have looked more closely at the welfare impact of the feedbackloop. This line of works points out to the fact that fire-selling assets during a crisis entails a negative externality, since the implied fall in asset prices tends to tighten the credit constraints of all the other agents -thereby making them more likely to also sell their own assets. As a consequence, a benevolent policymaker may be willing to restrict agents' borrowing ex ante in order to limit the risk of fire sale. This theory provides one possible justification for imposing a "macroprudential" regulation, in addition to the more traditional banking regulation.

General-equilibrium dampening of aggregate shocks

The discussion above illustrates the fact that general-equilibrium feedbacks can dramatically amplify the impact of "small" aggregate shocks, and stresses that such feedbacks can only be measured in the data by means of a theoretical model that incorporates them in the first place. But it is worth stressing that general-equilibrium effects do not necessarily manifest themselves through amplification: they can (2018). [START_REF]Frictional coordination[END_REF] wonders what classes of models can rationalize the Keynesian narrative that low aggregate demand can depress output. To frame the discussion in modern language, suppose that individual consumers suddenly value current consumption less than future consumption (say, marginal utility falls exogenously relative to future marginal utility). Holding prices constant, and aggregating over all consumers, this preference shock must translate into lower aggregate consumption demand, hence it is indeed a "negative demand shock". In a partial equilibrium setting with constant prices, this would translate into lower consumption. But prices cannot be considered constant in general equilibrium. As stressed by Angeletos ( 2018), in a Real Business Cycle model the drop in consumption generates an equal rise in savings and a fall in the real interest rate that boosts investment demand. Now let's take this reasoning one step further and assume that output uses labor only, so that there is no demand for capital on the part of firms. Still, if markets are complete households can potentially trade bonds between themselves, instead of lending capital to firms. But someone needs issue the bonds that the savers are willing to purchase, and no household hit by a negative marginal utility shock is willing to borrow. In general equilibrium, the (shadow) price of bonds must fall until households are again happy consuming the very same level of consumption that was planned before the marginal utility shock occurred: general-equilibrium adjustments in relative prices have completely eliminated the partial-equilibrium effect of the consumption shock. Of course, aggregate demand and output may fall after a consumption shock if prices are sticky and output is demand-determined; in this case actual output may fall below natural output, which is formally equivalent to a rise in monopolistic distortions [START_REF] Woodford | Interest and Prices[END_REF]. But here again, endogenous price adjustments, even muted, tend to reduce the direct (partial-equilibrium) effect of the consumption shock, except in the extreme case of constant nominal prices.

To summarize, partial-equilibrium intuitions or empirical evidence are uninformative about the likely effects of aggregate shocks, since we should expect those shocks to be amplified (due to strategic complementarities and feedback effects) or buffered (due to endogenous price adjustments) or both at once.

What Lessons for Macroeconomic Policy? The Example of Fiscal Multipliers

The importance of general-equilibrium effects implies that the impact of alternative macroeconomic policies cannot generally be estimated simply by extrapolating measures, however precise, based on "small", local policy shocks. The recent debates on the size of the fiscal multipliers, and notably the government spending multiplier, illustrate this point and deserve further discussion.

Formally, the government spending multiplier is defined as the growth in GDP induced by an exogenous increase in government spending scaled by GDP before the policy change. The empirical literature on this multiplier is considerable. Its main challenge is to measure the causal effect that goes from public expenditure to output, while many other mechanisms may affect the empirical correlation between these two variables. To make this point clear, imagine that government spending has no causal effect on output whatsoever. However, government spending varies systematically with output since it is higher in recession than in expansion (due to the automatic stabilizers), hence there is a reverse causality going from output to government spending. The endogenous response of government spending to output induces a negative correlation between these variables that can be wrongly attributed to a causal effect running from expenditure to ouput. In practice causality goes both ways, and moreover a number of confounding variables may correlate government spending and output independently of any causal link. In this context, how can one isolate the variations in government spending that are truly exogenous, in order to measure their causal effect on output?

The recent empirical literature has mostly relied on two distinct identification strategies to answer this question. The first strategy is to focus on a particular type of government spending shocks that are arguably not themselves caused by changes in GDP. The most common way to proceed is to consider as exogenous shocks the increase in military expenditure due to sudden, unanticipated deteriorations of the geopolitical context. These events generate variations in public spending that do not depend on the business cycle (although the cycle depends on it) and thus constitute in principle a valid basis for measuring the public spending multiplier. The multipliers obtained using this method vary between 0.8 and 1.2 for the United States (Hall, 2009[START_REF] Ramey | Macroeconomic shocks and their propagation[END_REF].

A second approach to estimating the causal effect from public spending to output relies on variations in local public expenditures to average out their nation-wide component. The study by [START_REF] Serrato | Estimating local fiscal multipliers[END_REF] provides a particularly telling illustration of this approach. Every ten years, the population of the United States is counted, so that the recorded population of each county changes. Following this, the federal government adjusts its financial allocation to adjust to these demographic changes: the counties whose population is revised upwards see their endowment increase, and those whose population is re-evaluated downward see it decrease. Unsurprisingly, this reallocation of federal funds between counties gives rise in each county to a variation in local public expenditure. By construction, these local variations are orthogonal to variations in US GDP as a whole as well as to other economy-wide factors (e.g., monetary policy) that are systematically related to GDP. Variations in local output generated by local variations in public spending thus form a valid basis for computing local spending multipliers. The authors find local multipliers close to 2, thus significantly higher than those obtained using macroeconomic data. The other studies adopting a similar approach also find high values of the multiplier, around 1.5. 4 In summary, the empirical literature on the public spending multiplier gives (for the United States):

-"nationals" multipliers between 0.8 and 1.2; -"local" multipliers between 1.5 and 2.

From a strictly empirical point of view, the estimation of local multipliers has two advantages over that of national multipliers. First, the exploitation of geographical disparities in public expenditure eliminates by construction any effect of the aggregate business cycle on public expenditure, which in principle offers a more reliable identification strategy than those based exclusively on macroeconomic data. Second, local multipliers tend to be more precisely estimated, partly because they make use of a much larger set of data.

One question that naturally arises is the reliability of these multipliers in terms of macroeconomic policy. Indeed, by their very nature, the spending shocks under study are of low amplitude, and partly offset each other from one region to another. Therefore, it is unlikely that these shocks will trigger the potentially powerful general-equilibrium effects of a large-scale shock at the level of a country as a whole. These general-equilibrium effects of public spending shocks have ambiguous consequences on the size of the multiplier: they can either lessen the direct microeconomic effects of the shock (for example, if the public expenditure shock is associated with a rise in the real interest rate, which reduces private expenditure), or instead amplify them (if for example one of the amplification mechanisms described in the previous section are set in motion). For the reasons explained above, we should expect both amplification and dampening effects to be simultaneously at work, leading them to partly offset each other. Unfortunately, this implies that local spending multipliers are of little help, on their own, in assessing the likely effect of a macroeconomic stimulus package -the ultimate question of interest. Does this mean that local multipliers are of no interest to macroeconomists? Not so. As shown by [START_REF] Nakamura | Fiscal stimulus in a monetary union: Evidence from US regions[END_REF] even if they do not directly inform us about the size of the aggregate multiplier, local multipliers (which they more accurately refer to as the "openeconomy relative multipliers") offer a powerful way of evaluating alternative macroeconomic models. For example, under the New Keynesian model, the size of the economy-wide public spending multiplier is conditional (like any fiscal policy) on the response of monetary policy: a strict targeting of inflation can lead the central bank to raise the real interest rate (via an increase in the path of policy rates) following a government spending shock, with the effect of reducing the observed fiscal multiplier; in contrast, a more accommodative monetary policy response would reinforce the expansionary impact of the fiscal stimulus. In as much as local multipliers are independent of economy-wide monetary policy, the alternative models make unambiguous predictions as to the size of their relative open-economy multipliers, which can then be compared to their empirically estimated counterparts. Nakamura and Steinsson show that this exercise leads to the rejection of the neoclassical Real Business Cycle model in favor of the New Keynesian model, which turns out to imply a much larger closed-economy aggregate multiplier. Their analysis makes it clear that it needs a fully specified macroeconomic model to turn estimated local multipliers into a policyrelevant, economy-wide government spending multiplier.

Concluding Remarks

In a letter to Harrod in 1938, in response to his presidential address to the Royal Economic Society, Keynes discusses the nature of economics in these terms: "It seems to me that economics is a branch of logic, a way of thinking; and that you do not repel sufficiently firmly attempts à la Schulz to turn it into a pseudo-natural-science. [...] Economics is a science of thinking in terms of models joined with the art of choosing models which are relevant to the contemporary world" (J.M. Keynes, Letter to Harrod, July 4, 1938).

We cannot better summarize what remains an essential characteristic of the study of business cycles and crises, namely the primacy of economic theory over empirical analysis. This remains true today even though the relationship between theory and data (and, more recently, micro data) is much tighter than when Keynes wrote these lines. This primacy of the theory makes the discipline necessarily more conjectural than other fields of economics, because statistical inferences are always conditional on relatively complex general-equilibrium models whose relative performance is difficult to evaluate. This is not a lack of scientificity, but rather the way in which the scientific approach manifest itself in this field of investigation.

Chart 1 .

 1 The liquidity trap and the deflationary spiral 2. See, for exemple,[START_REF] Eggertsson | The paradox of toil[END_REF] Krugman (2012), Christiano et al. (2015), andGust et al. (2015).

  introduced. Ravn and Sterk (

  dampen the impact of aggregate shocks, relative to what a naïve partial-equilibrium analysis (holding prices constant) could suggest. To illustrate this, let us push the line of argument of Angeletos

[START_REF] Eggertsson | The paradox of toil[END_REF] introduced the paradox of toil by studying the impact of labor supply shocks on equilibrium employment in a liquidity trap (he showed that a positive labor supply shock could actually lower employment, due to the inflationary impact of the shock on nominal wages and prices). Since then the same expression has been used to qualify the paradoxical effect of any supply shock on output in a liquidity trap.

See, for example,[START_REF] Acconcia | Mafia and public spending: Evidence on the fiscal multiplier from a quasi-natural experiment[END_REF] and[START_REF] Nakamura | Fiscal stimulus in a monetary union: Evidence from US regions[END_REF], as well as Fuchs-Schuendeln and Hassand (2016) for a survey.


F Plouraboué

H.-C Chang

Symmetry Breaking and Electrostatic Attraction between Two Identical Surfaces

By allowing the surface charge of one surface to affect the adsorption equilibrium of the other, we establish the existence of a long-range attractive interaction between two identical surfaces in an electrolyte containing polyvalent counter ions with a mean-field Poisson-Boltzmann approach.

A Stern electrostatic condition from linearization of the mass-action adsorption isotherm is used to capture how polyvalent ion condensation affects and reverses the surface charge. We furthermore establish a direct mapping between this Stern layer conditions and previously derived modified Mean-field formulations associated with correlated fluctuations theory. For a sufficiently potential-sensitive isotherm, anti-symmetric charge inversion can occur to produce an attractive force that increases with decreasing ionic strengths. Analyses of a mass-action isotherm produce force-separation relations, including an exponential far-field force decay distinct but consistent with previously proposed correlated fluctuation theories, and in quantitative agreement with experimental data.

I. INTRODUCTION

Recent atomic force microscopy studies have shown that, in the presence of a polyvalent counter-ion, two similarly charged or identical surfaces can develop an attractive force at a distance comparable to the Debye screening length λ [1,2]. This observation is most likely related to earlier reports on attraction between identical colloids in an electrolyte, although the role of poly-valency is not as well established for colloids [3,4]. It is also related to the condensation of likecharged molecules like DNAs [5]. It has been rigorously

shown that, for two identical spheres with constant potential [6] or constant surface field [7], the inter-particle interaction is necessarily repulsive according to the classical Poisson-Boltzmann (PB) mean-field theory [6]. For this reason, theories for like-charge attraction phenomena have sought mechanisms beyond the classical mean-field description to include spatial correlation of charge fluctuations [8][9][10][11].

Such fluctuation theories suggest cross-surface spatial correlation either between fluctuations of the condensed counter-ions in the Stern layer [12] or between surface Wigner crystals [13] that are formed by the same ions. However, when the two surfaces are separated by a distance larger than the Bjerrum length, thermal noise is expected to disrupt the fluctuation correlation of condensed ions on the two surfaces. Also strong concentration gradients develop within the Debye layer which may invalidate the assumptions of many correlated fluctuation theories.

Instead, one expects mean-field theory to be valid at separations comparable to the Debye length, when it is much larger than the Bjerrum length. Spatial correlation of fluctuations between the condensed charge and the surface charge on the same surface, which is over a distance comparable to or shorter than the Bjerrum length, affects the local screening and charge inversion [14]. This is particularly true for polyvalent counter-ions, which can cause local charge inversion and whose fluctuations are known to be important. Indeed, the twofluid model for free/condensed counterions proposed in [10] shows, with a modified mean-field theory, that the surface potential and the charge might be coupled by fluctuations on the condensed plates. Nevertheless this coupling was not previously considered in [10] where a zero potential is prescribed at the boundaries. Hence, long-range attraction over large Debye layers is quite distinct from bridging of both surfaces by the polyvalent counterion over short Stern layer separations [15].

Although AFM measurements of [1] and recent MD simulation of DNA condensation does not show surface ion ordering [16], correlated fluctuations should nevertheless play a role during ion condensation onto charged surfaces [17] especially for polyvalent electrolytes counterions [2]. However, connecting the fluctuations at a Bjerrum length separation from the surface to longer range attraction over the Debye length is still unclear. More specifically, theories for correlated layers charges between two parallel plates separated by a distance 2h much larger than the Debye length λ, give rises to an effective "long-range" dipolar attractive pressure, p ∼ 1/(2h) 3 [18], an almost dipolar one p ∼ ln[2h/λ]/(2h) 3 [8,10], or a modified exponential attraction p ∼ exp -4h/λ /2h [19,20]-also found in colloids interactions [21]-whereas the classical DLVO theory only predicts a repulsive exponential behaviour p ∼ exp -2h/λ . In this paper, we lump short-range correlated fluctuation effects into a linear empirical isotherm for polyvalent counter-ion condensation that captures charge inversion.

We then show that this linear Stern isotherm condition can induce asymmetric charge inversion and PB ion distributions due to nonlinear field screening/enhancement of the space charge and field-enhanced condensation effects on the two surfaces.

Moreover, we are able to show that, despite the simplicity of the model, the resulting force-separation relations display a variety of attractive/repulsive patterns which are quantitatively consistent with AFM measurements in the same polyvalent electrolyte. Fundamental to the theory is the observation that, for separation less than or comparable to λ, the field from one surface can affect the adsorption equilibrium and the Stern layer potential on the other. The field transmission across the large Debye layer separation with large concentration gradients must be described by a mean-field theory. For nearly neutral surfaces near the isoelectric point, where charge inversion is possible, it is shown that these two surfaces necessarily undergo opposite charge reversal to produce an attractive interaction. This attraction can be overwhelmed if the native charge density is too high to allow charge reversal on one of the surfaces. The paper is organized as follows. We describe in section II A the classical Stern-layer boundary condition in its usual formulation. We then establish in section II B a one-to-one mapping between Stern-layer parameters and a statistical mechanics "two-fluid" model previously proposed for fluctuation theory. In section III we first discuss the linear Debye-Hückel (D-H) approximation with Stern layers in order to find the bifurcation points for the non-linear Poisson-Boltzmann solutions which are further analyzed analytically in III B. Section IV is devoted to the numerical computation of the main caracteristics of these solutions such as the surface potential and the interaction force with Stern layer parameters and dimensionless gap. Finally our findings are compared with available experimental data in section VI. We first show in VI A that adsorption isotherm linearization leads to Stern layer parameter consistent with those obtained in section II B.

We then compare long-range and intermediate range force behaviour predicted by our theory with experimantal measurements in section VI B.

II. MODIFIED MEAN-FIELD DESCRIPTION

In this section we first described the classical Stern layer condition to account for ion condensation at the solid surface and its relation to adsoption isotherms. We then present in detail the statistical mechanical two-fluid model of Lau and Pincus [10] for counterion condensation. We then establish its equivalence, at the mean-field level, with a Sternlayer boundary condition, when considering Neumann Boundary condition rather than the previously considered homogeneous Dirichlet.

A. Stern layer boundary condition

Assuming rapid polyvalent counterion condensation kinetics, the isotherm stipulates that the surface charge density σ on a surface is a function of the Stern layer potential φ and the bulk electrolyte concentrations. However, since σ is proportional to the surface electric field, we obtain after linearization in φ at a particular reference state, the normal field

∂ z φ(±1) = ±K[φ(±1) -φ 0 ], (1)
where the inverse Stern slip length K is determined by various association/dissociation equilibrium constants of participating ions and φ 0 represents the compensate Zeta potential on the true surface due to polyvalent ion condensation. This value is exactly zero at the isoelectric point. In Stern's classical context, this boundary condition represents the field across the Stern layer with the condensed counterion and φ 0 is the potential just outside the Stern layer. The inverse Stern slip length K is then related to the thickness of the condensed layer and the surface density of the condensed counterions. The parameter K can take any sign, depending on the relative affinity of the surface to polyvalent and monovalent counterions. Our analysis is for two planar surfaces although extension to spheres can be carried out with the classical Derjaguin formulation. We scale the normal coordinate by h, the halfseparation between the two surfaces and the potential by RT /ZF to obtain an expression identical to (1) but with dimensionless variables. If we define the dimensionless separation squared as β = (h/λ) 2 , then the dimensionless K, scaled by the separation h can be decomposed in K = -µ √ β such that µ is only a function of the electrolyte composition.

The magnitude of µ measures the effect of the Stern layer potential (and external field) on the adsorption isotherm.

With the same scaling, the PB equation for symmetric electrolytes is

∂ 2 z φ = β sinh φ, (2)
which will be solved with (1) for z = ±1. We note that, (1) and ( 2) are symmetric to the transformation K → -K and z → -z. Hence the sign of K is arbitrary and will be taken to be negative, so that µ is a positive constant. For any value of φ 0 , symmetric solutions always exit, but anti-symmetric solutions are only permitted for φ 0 = 0.

B. Equivalence with "two-fluid" counterions condensation models

Let us first give here a brief description of the two-fluid model analyzed in [10,11,22].

The counterions of a charged solution are divided into a condensed fraction of charges of surface density n c and a prescribed surface charge density n 0 . The condensed charges are trapped onto the surface while the free charges can move into the bulk fluid. The surface density of the condensed counterions n c is not a priori prescribed but results from the balance between electrostatic interactions and fluctuations. In the frame of a point-like particles description on the surface the condensed counterions n c has been computed in [10,22] from the functional minimisation of the free energy which includes 2D and 3D entropic and electrostatic contributions. It has been found in [22] that correlation effects are responsible for a first order phase transition of the condensed phase, when varying the surface charge order parameter g = πZ 3 ℓ 2 B n 0 , where Z is the ion valency, ℓ B is the Bjerrum length defined by ℓ B ≡ e 2 /ǫk B T associated with the electron charge e, permittivity ǫ, Boltzmann constant k B and the temperature T . From the cubic dependence of g on Z, it is clear that the condensation is very sensitive to the ion valency, as quantitatively observed in many experiments. The minimisation of the free energy of the 2D and 3D entropic and electrostatic contributions leads to a modified mean field equation for the mean-field potential φ. This modified mean-field equation has been written for single counterion in equation ( 17) of [10] and straightforwardly extended to electrolyte solution in [11]. We write their modified mean field equation in dimensionless parameters for condensed counterions in a electrolyte placed between two plates located at position z = ±1:

∂ 2 z φ -β sinh φ = n R ℓ B h Z ± δ(z = ±1)+ πZ 2 n c ℓ B h ± δ(z = ±1)φ(z), (3)
where δ stands for the Dirac distribution. The two terms of the right-hand-side of ( 3)

have two distinct origins. The first one is associated with the point-like contribution of the condensed charges onto the surfaces which naturally result from an electrostatic contribution.

It is nevertheless interesting to note that the prefactor of this term is not the obvious imposed surface charge density n 0 but the bare charge n R which results from the renormalyzed bare surface charge density n R = (n 0 -Zn c ). The second term result from the fluctuation contribution to the condensate and couples the mean-field potential to the condensate surface charge. The contribution of this second term to the surface interaction has been previously

ignored in [10,11,22] from imposing a zero potential boundary condition φ(±1) = 0 for the PB mean-field φ. Hence, from using a Dirichlet boundary condition which annihilates the contribution of the second term of Eq. ( 3) right-hand-side previous studies have then prescribed a zero contribution of the fluctuation dependence of the mean-field solutions.

They subsequently analyzed the challenging estimate of fluctuations contribution to the attraction between surfaces, whilst discarding their influence at the mean-field level.

Our contribution in this paper is to show that keeping with the modified mean-field

description
-∂ z φ, ∂ z Φ -β sinh φ, Φ = n R ℓ B h Z [Φ(-1) + Φ(1)]+ πZ 2 ℓ B n c h[φ(-1)Φ(-1) + φ(1)Φ(1)] (4)
It can be seen in the weak formulation that the source terms are essentially associated with potential value at the surface φ(±1). Now, using the same weak formulation for PB equations (2) but prescribing instead the Stern layer boundary conditions (1) leads to the following problem:

-∂ z φ, ∂ z Φ -β sinh φ, Φ = µ √ βφ 0 [Φ(-1) + Φ(1)]+ µ √ β[φ(-1)Φ(-1) + φ(1)Φ(1)] (5)
Hence, formulations (4) and ( 5) are formally equivalent. They can be exactly mapped one onto another provided that the parameters µ and φ 0 fulfill

µ = πZ 2 ℓ B λn c (6)
φ 0 = n R n c πZ 3 (7)
Relations ( 6) and ( 7) then provides a statistical mechanics based Stern layer formulation that we have not been able to find elsewhere in the previous literature. Parameters µ and φ 0 then depends on the Bjerrum length ℓ B , the ions valency Z, the condensed surface charge n c and the effective bare surface charge n R . n c and n R also depends on the order parameter g = πZ 3 ℓ 2 B n 0 , through a first order transition and have been computed numerically and asymptotically in certain limits in [22].

In section VI A we will nevertheless prefer to derive these parameters from the adsorption isotherm linearization parameters, which are more easy to relate to some given experimental condition.

III. SOLUTION TO THE PB PROBLEM WITH STERN LAYERS

In this section we solve the PB problem associated with boundary conditions (1). We first investigate the linear approximation associated with the Debye-Hückel (D-H) approximation, and then explore the features of the solution to the full non-linear problem.

A. D-H approximation

Even though the linearized version of (2) does not produce any interaction, we still use the linear Debye-Hückel (D-H) approximation to determine the bifurcation of non-trivial solutions. Two linearized solutions exist, one symmetric (denoted with a capital S) and one non-symmetric (denoted with NS) with respect to reflection across midplane,

φ S = cosh( βz), (8)
φ N S = e √ βz 2 + ( 2µφ 0 µ 2 -1 -1) e - √ βz 2 . ( 9
)
The coefficients of these eigenfunctions cannot be specified for the homogenized versions of ( 1) and ( 2). However, the corresponding linear operator becomes singular at

β S c = ln -µφ 0 + µ 2 φ 2 0 + (µ 2 -1) (1 -µ) 2 , ( 10
)
β N S c = 1 2 ln µ + 1 µ -1 2 , (11)
where non-trivial solutions can bifurcate from the trivial solution 0. We note that for nonsymmetric solutions (anti-symmetric in the case φ 0 = 0), only |µ| > 1 is possible, or, more precisely, µ > 1 since we are restricted to positive µ values.

B. Solution to Poisson-Boltzmann problem

Equation ( 2) is integrable and its first integral is

1 2 (∂ z φ) 2 = β cosh φ + d, (12)
where d is a constant which depends on µ, β, φ 0 and the surface potential. Evaluating this integral at z = ±1 as a reference point, one obtains

d ′ = d β = 1 2 [µ(φ(±1) -φ 0 )] 2 -cosh φ(±1). (13)
This parameter can be related to the force between the two surfaces. Using Green's theorem and including the osmotic pressure from the bulk solution at the infinities behind the two surfaces [6], the pressure (force per unit area) between the two surfaces is can be derived, where the dimensionless pressure p has been normalized by the Debye Maxwell stress ǫ(RT /zF λ) 2 for convenience to remove the dependence on separation. It is independent of z and hence its value at z = 1 is evaluated by using (13). Looking at the last expression of (14) it is obvious that the osmotic pressure difference in the second term is always repulsive and the Maxwell pressure in the first term is always attractive, leading apparently to possible attractive or repulsive forces depending on φ(1) or d ′ . Nevertheless it can be easily shown by evaluating the pressure at z = 0 for the case of φ(0) = 0 associated with anti-symmetric solution for φ 0 = 0, in the first expression of ( 14) that the only remaining term at the origin is the attractive Maxwell pressure. Hence all anti-symmetric solutions near the isoelectric point are attractive. Furthermore, some observations can be made at the (D-H) limit of small φ(1). Expansion of the osmotic pressure term shows that attractive interaction is only possible for µ > 1, which is consistent with the previous observation that anti-symmetric solutions are only possible for this parameter range. The required solution can be simplified by the Boltzmann transformation

p = -1 2β (∂ z φ) 2 + cosh φ -1 = -d ′ -1 p = -[µ(φ(1)-φ 0 )] 2 2 + cosh φ(1) -1, (14)
ψ = exp -φ 2 , (15)
such that the first integral ( 12) is transformed to

∂ z ψ = ± √ β 2 ψ 4 + 2 d β ψ 2 + 1 ∂ z ψ = ± √ β 2 (ψ 2 -α -)(ψ 2 -α + ), (16)
where,

α ± = -d ′ ± d ′ 2 -1, (17)
are complex conjugate roots. Integrating ( 16) one more time leads to

± √ βz 2 + c = - √ α -F 1 (ψ √ α + , α -) , (18)
where F 1 is the elliptic integral of the first kind. Evaluating (18) at z = ±1 gives a transcendental equation for the potential values at boundaries:

β = - √ α -[F 1 (ψ(1) √ α + , α -) -F 1 (ψ(-1) √ α + , α -)] , (19)
which complements Eq.( 13), (15) and Eq.( 17) to produce a set of two transcendental equations for φ(±1) that can only be solved numerically.

IV. NUMERICAL COMPUTATION AND RESULTS

Albeit α ± is complex, and thus elliptic integral of the first kind with complex variables are needed for this computation, the solutions are always real and are symmetric with respect to α ± . We use a Newton-Picard method for the evaluation of the unknown surface field φ(1), and continuously vary the parameters β and µ with a continuation method. The pitchfork bifurcation of these antisymmetric solutions from the bifurcation point β N S c (µ) of ( 11) is shown in Fig. 1 for different values of µ. It is evident that φ(1) first rises beyond a critical dimensionless separation β N S c and very rapidly reaches a constant asymptotic value φ ∞ (µ) which depends only on µ at large separations, i.e for large value of β. This asymptotic value is concomitantly reached when d ′ tends to -1, and α ± to 1 so that the arguments of the elliptic function in ( 18) reaches 1 where it display a unique logarithmic singularity to compensate for large value of β on the left-hand-side. In this limit, one can evaluate the potential φ ∞ (µ) from solving the simple transcendental equation

[µ(φ ∞ -φ 0 )] 2 = 2(cosh φ ∞ -1) (20)
found from ( 13). This values is illustrated for φ 0 = 0 in Fig. 1. Furthermore, from knowing this value, it is easy to find the anti-symmetric solution which is depicted for half of the domain in Fig. 2 for µ = 2.

One can see from this plot the unconfined charge distribution and long-range attraction of anti-symmetric solutions. For small values of β, the potential gradient at the origin z = 0 monotonically increases with β. For a critical ratio of dimensionless half-gap β m (e.g equal to β m ≃ 0.391 when µ = 2), this gradient at the origin reaches a maximum, after which the potential profile starts to deform non-linearly. Increasing β further more both decreases the potential gradient at the origin and increases it at the wall. When the surface potential at z = 1 reaches the asymptotic value φ ∞ (µ) solution of (20), the potential gradient at the wall monotonically increases with √ β according to (1).

A boundary layer develops near the wall with large potential gradient whilst the potential rapidly decreases to zero away from the wall. The surface charge, which is proportional to the potential gradient, is then increasing as √ β leading to an asymptotic attractive force/potential at long range that will be subsequently studied in next section.

V. ATTRACTIVE BEHAVIOUR

From knowing the surface potential one can easily deduce from ( 14) the negative pressure associated with iso-electric point anti-symmetric solutions. This negative pressure is plotted in Fig. 3a. The sharp maximum of the pressure in Fig. 3 the parameter µ. When moving further apart from the iso-electric point, one find that the attractive region exists only for a finite range of separation β which are larger than β N S c , as illustrated in Fig. 3b. Repulsion due to native charges responsible for the effective Zeta potential φ 0 reduces the attractive region until it disappears completely beyond a critical φ 0 . The repulsion also gives rise to a shallow maximum, a threshold pressure, as seen in Fig. 3b. This leads to the "phase-diagram" in Fig. 4a for the attraction/repulsion on the β -φ 0 plane. Attractive interaction occurs near isoelectric points with minimum surface charge (low φ 0 ) and with high sensitivity of the condensation isotherm to external field (high µ).

When attractive interaction is possible, a shallow threshold barrier exists except for nearly isoelectric conditions.

As previously discussed, the large gap limit is interesting to consider and compare to the usual (D-H) limit for the force. In this limit, the parameter d ′ approaches value -1, which corresponds to the vicinity of the logarithmic singularity of the elliptic function. Using known asymptotic behavior of the elliptic function near this singularity, it is possible to evaluate analytically the force distribution using, ( 13), ( 14) and ( 18)

for β ≫ 1 p ∼ -exp -2 √ β ∼ -exp -2 h λ (21)
We compare this asymptotic behavior against the numerical results in the inset of Fig. (3)a and show very good agreement for β values as small as 3. This asymptotic force is also consistent with our numerical values for both attractive or repulsive long-range forces when φ 0 = 0 (not shown). It is interesting to note that the slope of the curves is independent of µ, as predicted by ( 21) and is consistent with the constant asymptote of the surface potential at large separation in Fig. 1.

VI. COMPARISON WITH EXPERIMENTS

A. Derivation of Stern layer parameters from adsorption isotherm linearization

In this section we use the "classical" derivation for Stern layer parameters from adsorption isotherm linearization. We evaluate the µ and φ 0 parameters from the isotherm proposed in Zohar et al. [1] for the net surface charge density σ at different trivalent Cobalt hexamine concentrations. The isotherm is derived from simple mass-action kinetics whose parameters are either measured or estimated with simple arguments, shown to be inconsistent with phase-locked Wigner crystals

σ(φ) = -n Si e 1 -K 1 [Co] exp -3φ 1 + K 2 [Co] exp -3φ (22)
Where e is the electron charge, [Co] the Cobalt hexamine concentration, φ the dimensionless potential rescaled by RT /ZF , K 1 = 10 4 M -1 and K 2 = 5 10 3 M -1 are two empirical constants related to the association and dissociation rates and n Si is the density of ionizable silanol groups which is estimated to range between n Si ≃ 0.3-1 nm -2 . The surface electro-neutral condition σ(φ 0 ) = 0 leads to a simple dependence between the compensated Zeta potential φ 0 at a particular counterion concentration

φ 0 = ln (K 1 [Co]) 3 , (23)
A first order Taylor expansion of ( 22) about this reference point leads to

σ(φ) ≃ -2n Si e(φ -φ 0 ) (24)
The application of Gauss theorem at the surface location z = 1, for dimensionless potential rescaled by RT /ZF = k B T /Ze leads to

∂ z φ = -2 Ze k B T n Si eh ǫ (φ -φ 0 ) (25)
From identifying the constant of this affine relation to the Stern condition (1) with K =

-µh/λ one finds

µ = 2Zℓ B λn Si (26)
where, again, ℓ B ≡ e 2 /ǫk B T is the Bjerrum length. It is then interesting to compare (26)

with the previous expression derived from the theoretical two-fluid model of Lau et al. in (6). The condensed ion density n C in the theoretical derivation of constant µ proposed in ( 6) is replaced by the density of ionizable silanol groups n Si in its empirical counterpart given in (26). Since they both represent the surface density of negative charges at the solid surface, this give a consistent picture of how the parameter µ is related to physical constants.

Moreover, since n C is a theoretical prediction of the surface charge density (which could be related to other physical parameters from charge condensation mecanism as explained in [22]) where n Si is related to an experimental estimate of the same quantity, both expression are fully consistent with either a theoretically or experimentally derivated expression for µ. Finally, it is hence surprising that such different approaches for estimating parameter µ only quantitatively differ by a prefactor which is 2Z in one case and πZ 2 in the other case. The different dependence with the valency Z could be attributed to the fact that in the experimental case, the negative charges associated with silanol groups are generally considered as monovalent.

The range of possible experimental value for n Si , leads to µ = 1.32-4.4. Furthermore relation (23) produces negative values of φ 0 for these experimental condition. However, as the problem is invariant with the transformation (z, φ, φ 0 ) → (-z, -φ, -φ 0 ), our results obtained for positive φ 0 also applies for negative ones.

B. Comparison with experimental force

The configuration which is analyzed in this paper is associated with two parallel identical surfaces. Nevertheless, the available experimental fine measurements of the force between similar surfaces in an electrolyte have all been performed with Atomic Force Microscopy (AFM) [1,2,4]. Those measures are associated with a sphere/plane configurations. Nevertheless some features of the predicted attraction force between two parallel planes can be possibly extended when the ratio between the AFM tip distance to the wall to its radius of curvature is small [4,20,23]. A first possible direct comparison between theory and experiments can be found from analysing the far-field behavior of the force which in both cases should be the same up to some irrelevant prefactors [23].

Long-range behaviour

As mentionned earlier in section I, different theoretical predictions have been proposed for the long-range behaviour of the possibly attractive force of an electrolyte confined between two identical walls distant by a half-gap h. Nevertheless, few of these predictions have been sucessfully compared with experiments but the modified exponential attraction

p ∼ exp -4h/λ 2h ( 27
)
has been proposed in [19,20] and experimentally tested in [4].

Nevertheless, there might be some caution for a definite answer to this question since the proposed modifyed exponential behavior very slightly differ from a pure exponential behaviour more recently reported in [1,2]. Nevertheless, to our knowledge, the possibility of non exponential, algebraic far field decay proposed in [8,10,18] has never received any convincing experimental evidence. Futhermore, since most of the experimental evidence for a modified exponential attraction is tested on semi-logarithmic coordinates in [4], and given the fact that experimental observation are also relatively noisy, it is hard to draw a definite conclusion on a possible logarithmic corrections to a mere linear trend of ln p versus h in (27).

Moreover, most of the correction to the linear behaviour in figure 2 of [4] is concentrated in the close region where it is higly possible that the far-field asymptotic regime break-down.

On the contrary, our exponential decay of long-range attraction ( 21) is consistent with various experimental observations [1,2]. To the best of our knowledge, our approach is the only one which predicts a pure exponential asymptote for the attractive force found in [2].

We verify this statement by analysing the slope of the exponential decay of the experimental measurements of [2]. To compare more quantitatively the current theory between two planar surfaces to that between an AFM tip (an attached spherical bead) and a flat surface in the experiment, an averaging over the AFM tip is carried out. The justification of this approximation in the framework of a rigourous asymptotic lubrication approach can be found in [23]. Furthermore, the computation of the macroscopic Force which is experienced by the AFM tip necessitates the knowledge of the local pressure dependence p(h(r)) versus the local distance h(r) between the plane surface and the AFM tip. In principle this dependence being non analytic, its computation necessitates some numerical integration as performed in [23].

A shortcut to this more complete computation is to find some simple expression for the pressure dependence with h(r) which reasonnably fits the far field behavior of the pressure dependence with h displayed in figures 3.

A very good estimate of the pressure dependence with dimensionless half-gap β is p(β) ≃ β -5/2 = (h/λ) -5 as illustrated in Fig. 3b.

Averaging this pressure from the true tip separation h m to infinity where the interaction vanishes and setting it equal to the pressure of an equivalent separation between two planes,

P ( h /λ) = h λ -5 = ∞ hm λ h -5 dh = 1 4 h m λ -4 (28)
one gets a dimensionless equivalent half-gap separation of

h λ = 4 1 5 [ h m λ ] 4 5 , (29)
whose variation with the initial minimum distance h m is almost linear. Using the equivalent gap h deduced from the observed h m from (29), associated with LaCl 3 salt analyzed in [2], we plot in figure 5 the experimental pressure as well as our theoretical prediction p ∼ exp -2 h /λ . The resulting slope in log-linear scale very satisfactorily compared with our theoretical prediction as previously stated in relation (1) of Besterman et al.'s paper [2] independently of parameter µ. Hence not only our prediction quantitatively agree with a pure exponential decay, but the obserbed slope in semi-logarithmic representation also leads to satisfactory comparison with experimental data.

Intermediate distance force

We can carry-on this approximation within intermediate separation distances as well as when varying the concentration conditions. To make a possible link between some experimental conditions and theoretical predictions, one needs to evaluate the Stern layer parameters µ and φ 0 . This is not possible with Besterman's data since no absorption isotherm is given for their experimental conditions [2] but this is possible with's Zohar et al.s data as was done in section VI A.

We choose the cross-section area of the spherical AFM tip at a separation of two Debye lengths as a reference area to calculate the force such that

F = 4πRλǫ k B T Zeλ 2 p ( 30
)
with the AFM tip radius R = 510 -6 m, the liquid permitivity ǫ = 80ǫ 0 , the temperature 293K, represented in Fig. 6. They display the same qualitative trends and variations as the measured ones in Fig. 2 of [1] -and are quantitatively consistent within the same order of magnitude. Hence, our predictions not only lead to an accurate estimate of the asymptotic long-range decay of the force, but it also gives consistent intermediate distance behaviours.

VII. CONCLUSION

We have shown that it is possible to obtain an attractive anti-symmetric solution for the potential distribution of an electrolyte between two identical surfaces by taking into account the coupling between surface charge across the double layer due to the effect of each field on the condensation of polyvalent counterions on the other, as captured by a simple linear Stern condition derived from common isotherms obeying mass-action kinetics. We show analytically that asymmetric charge inversion and attraction are possible and we produce numerical results that are qualitatively and quantitatively consistent with experimental data.

Our results reconcile correlated fluctuation theory based attraction, Stern layer models and available experimental observations.

1 - 1

 11 is enough to satisfactorily describe most of the experimental observations on charged like attraction. We now first show how the modified mean-field description (3) is formally equivalent with the Stern layer boundary condition (1). Instead of considering the strong formulation (3), let us now establish the weak formulation of this non-linear elliptic problem associated with the dual test functions Φ. Instead of considering an homogeneous Dirichlet boundary condition assiciated with (3), we look for solutions associated with an homogeneous Neumann boundary condition at the surface ∂ z φ(±1) = 0. Using notation φ, Φ ≡ φΦdz, and integrating the weak formulation of (3) by parts one gets:

FIG. 1 :

 1 FIG.1: Surface potential of asymmetric solution φ(1) at z = 1 versus dimensionless gap β for φ 0 = 0. The dotted lines display the value φ ∞ (1) obtained from solving(20).

FIG. 2 :

 2 FIG. 2: Potential field versus distance z for different dimensionless gap values, µ = 2 and φ 0 = 0.

FIG. 3 :

 3 FIG. 3: (a) Semi-logarithmic plot of the attractive force per unit surface p versus the square of the dimensionless half-gap β for different values of parameter µ and φ 0 = 0. Inset: the corresponding large-gap asymptotic prediction (21) is shown with straight-lines for each value of µ.(b) Same convention in usual coordinates for µ = 2 and three values of the compensate potential φ 0 . The continuous curve shows a -β -5/2 fit to one of the pressure variation.

FIG. 4 :

 4 FIG. 4: Investigation of attractiveness in the phase-space of the parameters of the problem. (a) (β -φ 0 ) plane for µ = 2. R stands for repulsion and A for Attraction. Point P (β ≃ 0.465, φ 0 ≃ 0.9) corresponds to the minimum equal to zero of the pressure curve represented on figure 3b (b) (µ-φ 0 ) plane. Long-dashed lines locate point P in this plane.

FIG. 5 :

 5 FIG.5: Test of the long-range scaling(21) for the data collected from Fig.1bof[2] associated with different concentrations of LaCl 3 salt with an adjusted equivalent half-mean separation h . The experimental Debye-length values λ are those given in[2].

FIG. 6 :

 6 FIG.6: Plot of the predicted dimensional force F (30) for λ = 10nm, three values of Cobalthexamine concentration relative to NaCl 1/200, 1/100, 1/25 which corresponds to 4.510 -6 M, 910 -6 M, 9.610 -5 M and gives the following compensated Zeta potential φ 0 = -1, -0.8, -0.013 analyzed in Fig.2of[1], using a density of ionizable silanol groups n Si = 0.48nm -2 leading to µ = 2.14. No adjustable parameter has been used.

Acknowledgments

We are grateful for inputs from Y.E. Zhu, X. Cheng, C. Beaume and S. Basuray, and fruitful discussions with Pr. E. Trizac.


PhD Jacques Dainat
email: jacques.dainat@gmail.com

P Pontarotti

Double space

Keywords: Pseudogenes, phylogeny, speciation, orthology, paralogy

The discovery that pseudogenes are involved in important biological processes has excited enthusiasm and increased the research interest on them. An accurate detection and analysis of pseudogenes can be achieved using comparative methods, but only the use of phylogenetic tools can provide accurate information about their birth, their evolution and their death, hence about the impact that they have on genes and genomes. Here, phylogenetic methods that allow studying pseudogene history are described.

Introduction

Genome analyses have shown the presence of numerous pseudogenes in all organisms. For example, Mus musculus and Homo sapiens both contain around 20 000 pseudogenes [START_REF] Bischof | Genome-wide identification of pseudogenes capable of disease-causing gene conversion[END_REF][START_REF] Khelifi | HOPPSIGEN: a database of human and mouse processed pseudogenes[END_REF][START_REF] Torrents | A genome-wide survey of human pseudogenes[END_REF][START_REF] Zhang | Comparative analysis of processed pseudogenes in the mouse and human genomes[END_REF]. Pseudogenes are categorized into three main groups according to the underlying mechanisms involved in their emergence: processed pseudogenes, unprocessed pseudogenes and unitary pseudogenes (Fig. 1).

• The first are retrotransposed pseudogenes. They are considered as dead-on-arrival [START_REF] Vanin | Processed pseudogenes: characteristics and evolution[END_REF], having suffered lethal damage during the copying process or missing important flanking regions [START_REF] Gerstein | The real life of pseudogenes[END_REF][START_REF] Satta | Primate Evolution: Gene Loss and Inactivation[END_REF][START_REF] Wang | Gene losses during human origins[END_REF]. • The second are duplicated pseudogenes. They originate from a duplication event of a functional gene. Gene duplication can represent an important source of diversity for organisms [START_REF] Fischer | Impact of recurrent gene duplication on adaptation of plant genomes[END_REF][START_REF] Magadum | Gene duplication as a major force in evolution[END_REF]. There are four potential outcomes for a duplicated gene: o Not subject to any selection pressure, it accumulates mutations and becomes a pseudogene. • The third group of pseudogenes is called unitary pseudogenes because they have no functional orthologous counterpart [START_REF] Mitchell | Inferring the pattern of spontaneous mutation from the pattern of substitution in unitary pseudogenes of Mycobacterium leprae and a comparison of mutation patterns among distantly related organisms[END_REF]. They originate from the pseudogenization of a gene that has been functional and subjected to selection for a long time period. Indeed, orthologs of such unitary pseudogenes still exists in other species. They are of special interest because their losses may be linked to functional losses.

Non-functional pseudogenes are not subjected to any selection pressure and drift until becoming unrecognizable [START_REF] Li | Pseudogenes as a paradigm of neutral evolution[END_REF]. They eventually become part of the genomic noise, formerly called "junk" DNA. The divergence rate of pseudogenes in avian and human lineages is estimated to be around 2% by million years [START_REF] Nachman | Estimate of the mutation rate per nucleotide in humans[END_REF][START_REF] Weir | Calibrating the avian molecular clock[END_REF], therefore the timeframe of a complete disappearance of a pseudogene can be estimated to 50 million years. If negative selection is in place, as it might be assumed in the case of the "less-is-more" hypothesis [START_REF] Olson | When less is more: Gene loss as an engine of evolutionary change[END_REF], the pseudogenization process can be even faster. Consequently, after a certain divergence time, it might become impossible to identify the pseudogene on the basis of similarity to its orthologs.

On the contrary, pseudogenes that maintain a transcriptional or even translational activity and are involved in biological processes [START_REF] Chan | Transcribed pseudogene ψPPM1K generates endogenous siRNA to suppress oncogenic cell growth in hepatocellular carcinoma[END_REF][START_REF] Hirotsune | An expressed pseudogene regulates the messenger-RNA stability of its homologous coding gene[END_REF][START_REF] Wen | Pseudogenes are not pseudo any more[END_REF][START_REF] Pink | Pseudogenes: pseudo-functional or key regulators in health and disease?[END_REF] can be preserved by positive selection and evolve like genes. Currently, the majority of methods to study the life of pseudogenes are based on comparative genomics [START_REF] Zhang | Identification and analysis of unitary pseudogenes: historic and contemporary gene losses in humans and other primates[END_REF][START_REF] Zhu | Comparative genomics search for losses of long-established genes on the human lineage[END_REF], which analyse sequence similarity and synteny (physical colocalization of multiple genetic loci on the same chromosome) to deduce the relationship occurring among different sequences. However, when the pseudogenization process has to be studied on a large evolution scale (superior to some million years), it should be considered that the structure of pseudogenic sequences can change compared to that of their parental genes, and the same may apply to synteny. These changes can in turn induce analysis difficulties and miss-interpretations and render comparative genomics a suboptimal approach.

During the last years, the increase in calculation capacity has enabled the use of phylogenetic methods to study robustly and rapidly variable biological entities such as pseudogenes. The phylogenetic methods are more advantageous than the classical methods of comparative genomics and particularly suitable to elucidate pseudogene birth, age and loss. This is because they allow inferring more robust relationships between the studied sequences. For example, a missing gene in a phylogeny is likely to be caused by a pseudogenization event [START_REF] Costello | Limitations of pseudogenes in identifying gene losses[END_REF]. By deciphering the phylogeny, it is possible to deduce in which ancestor the gene was present and starting from which ancestor the gene was lost. The analysis of these features (presence/absence) by the parsimony method [START_REF] Farris | Phylogenetic Analysis Under Dollo's Law[END_REF][START_REF] Mirkin | Algorithms for computing parsimonious evolutionary scenarios for genome evolution , the last universal common ancestor and dominance of horizontal gene transfer in the evolution of prokaryotes[END_REF][START_REF] Sankoff | Locating the vertices of a Steiner tree in an arbitrary metric space[END_REF][START_REF] Sankoff | Minimal mutation trees in sequences[END_REF], allows determining the appearance period of the pseudogene, which corresponds to the period of gene deactivation. Another case that exemplifies the use of phylogeny in the study of pseudogenes is represented by duplication events occurring before a speciation. Such events give rise to co-orthologs. The pseudogenization of one of the two duplicated genes should not lead to the loss of the function of the gene, due to the presence of a still functional copy. Consequently, this kind of clue allows excluding the hypothesis that the pseudogene is of the unitary kind.

In this chapter, we describe how phylogenetic methods can be applied to the identification and characterization of pseudogenes.

Methods

Pseudogene identification using phylogenetic tools

To identify pseudogenes the strategy is to look for sequences similar to functional genes and containing disrupting mutations. It can be performed "manually" using the sequence of the gene of interest as query and searching for similar sequences with a Blast tool (BlastP, TBlastN, BlastZ, etc) [START_REF] Altschul | Basic local alignment search tool[END_REF][START_REF] Schwartz | Human-mouse alignments with BLASTZ[END_REF]. Among the retrieved sequences, those not annotated as gene or/and having deleterious mutations that disrupt the open reading frame (ORF) can be considered as pseudogenes. Some pipelines automate that principle [START_REF] Khelifi | HOPPSIGEN: a database of human and mouse processed pseudogenes[END_REF][START_REF] Ortutay | PseudoGeneQuest -service for identification of different pseudogene types in the human genome[END_REF][START_REF] Zhang | PseudoPipe: an automated pseudogene identification pipeline[END_REF]. This approach does not allow to determine precisely the relationship between the sequence used as a query and the putative pseudogenic sequences detected. Excluded rare cases of convergence, those sequences are either orthologs or paralogs. This information is crucial for meaningful interpretation of the pseudogenization process. Moreover, there is no way to infer with precision temporality of the event.

A phylogenetic approach has the benefit to determine precisely the relationship of the studied sequences and the relative period when a pseudogenisation event has occurred. This relative period is determined by the speciation event of the species used for the phylogeny (see Notes 1).

The first step for a phylogenetic approach would naturally consist in building a phylogenetic tree of the gene or gene family of interest (see Notes 1-4).

Then, reading the phylogeny we can detect the missing genes and consequently the putative gene loss events (see Note 5). The putative losses are found by comparing the species tree of the chosen species with the phylogenetic tree of the gene of interest in those species (Fig. 2). The phylogenetic tree analysis can be done using software such as Phylopattern, Bio::Phylo and DendroPy [START_REF] Gouret | PhyloPattern: regular expressions to identify complex patterns in phylogenetic trees[END_REF][START_REF] Sukumaran | DendroPy: a Python library for phylogenetic computing[END_REF][START_REF] Vos | BIO::Phylo-phyloinformatic analysis using perl[END_REF].

When a gene is detected as missing in the genome of a species, the third step consists in searching for putative homologous sequences of that gene by a genome-wide BLAST. To optimize the chance to detect the sequence, the Blast should be performed using an orthologous gene retrieved from the phylogeny and such gene should be as close as possible to the investigated species. In the example reported in Figure 2, the loss of the gene of interest in H. sapiens should be investigated using the orthologous gene in M. mulatta. When no Blast hit is found, considering the genome assembly of the species as complete we can conclude that the gene has likely become genomic noise. But the available genome assemblies have variable numbers of missing regions, so any interpretation should be cautious. On the contrary, when a hit is found, its homology needs to be verified.

In the fourth step, we try to confirm the homology of the sequence retrieved in the previous step. The verification can be performed either by a synteny study or by a phylogeny study. To be orthologous by synteny, the retrieved sequence has to be located at the same position in the genome as the orthologous sequence use as query. Furthermore, the neighbouring genes should be the same. However, if distant species are compared in the synteny study (old speciation event), it is possible that the synteny has changed during evolution. In this respect, the phylogeny study ensures a higher grade of accuracy in inferring orthology, while excluding paralogy, or no homology (just similarity). This can be of special importance when looking at unitary pseudogenes (Fig. 3). To verify the relationship of the retrieved sequence with the sequence used as query, a new phylogeny must be built (see Note 6). Then, the orthologs of the retrieved sequence must be identified in this new phylogeny and compared with the orthologs defined in the old one (see Note 7). This can be done using software allowing reading phylogenies [START_REF] Gouret | PhyloPattern: regular expressions to identify complex patterns in phylogenetic trees[END_REF][START_REF] Sukumaran | DendroPy: a Python library for phylogenetic computing[END_REF][START_REF] Vos | BIO::Phylo-phyloinformatic analysis using perl[END_REF]. If both analyses give common orthologs, then the orthology of the Blast hit sequence is confirmed (see Notes 8 and 9).

Once the homology is defined, the Blast hit sequence should be subjected to in depth study (see Note 10). Two scenarios can occur. In the first case scenario, the sequence is already annotated within the genome: i) it is a gene present in the database used, but its sequence was not kept in the phylogeny building process; ii) it is a gene not present in the database used, but present in other databases; iii) it is indeed a pseudogene already described (see Note 11). In the second case scenario, the sequence is not annotated: i) it is a gene forgotten during the genome annotation process; ii) it is indeed an unannotated pseudogene. In reaching these conclusions, a more thorough investigation must be carried on (e.g., a gene prediction, ORF sanity, dN/dS ratio, use of RNAseq to check transcription and the sanity of the transcribed sequence). When a deleterious mutation is detected, it can be concluded that the sequence is a pseudogene. For ancient pseudogenes, the disruptions are easy detect, but this can be more difficult when only few disruptive mutations are presents. This can be particularly tricky because genes exist with stop-codon read-through [START_REF] Jungreis | Evolutionary dynamics of abundant stop codon readthrough[END_REF][START_REF] Stark | Discovery of functional elements in 12 Drosophila genomes using evolutionary signatures[END_REF][START_REF] Loughran | Evidence of efficient stop codon readthrough in four mammalian genes[END_REF]. The mutations searched for can be of various types: micro-mutations, such as substitutions or indels that engender the appearance of premature stop codons in the ORF (Fig. 4), but also macro-mutations that engender the loss of entire exons of the original sequence. Macro-mutations can be the result of a series of micro-mutations or of a recombination event (see Notes 12 and 13).

Pseudogene dating and characterization using phylogenetic tools

Once a pseudogene has been identified, the phylogeny approach allows getting insights into its evolutionary history and in particular on the time frame of the pseudogenization events. If there are several pseudogenes in the gene family studied, putting them together in the same phylogeny can improve the accuracy of the study. Thanks to the use of parsimony methods [START_REF] Farris | Phylogenetic Analysis Under Dollo's Law[END_REF][START_REF] Mirkin | Algorithms for computing parsimonious evolutionary scenarios for genome evolution , the last universal common ancestor and dominance of horizontal gene transfer in the evolution of prokaryotes[END_REF][START_REF] Sankoff | Locating the vertices of a Steiner tree in an arbitrary metric space[END_REF][START_REF] Sankoff | Minimal mutation trees in sequences[END_REF], it is in fact possible to reconstruct the ancestral state starting from the contemporary states that are present on the leaves of the phylogeny. In the example reported in Figure 5, the parsimony reconstruction allows to define three pseudogenization events for four pseudogenes. The first pseudogenization event (P1) occurred in R. norvegicus after the split with M. musculus and produced a unitary pseudogene that has no functional co-ortholog counterpart. The second pseudogenization event (P2) is specific of M. mulatta. The third pseudogenization event (P3) is shared between H. sapiens and M. mulatta and concerns the Catarrhini phylum. Further analyses can be performed to infer the retrotransposed or duplicated nature of the pseudogenes produced by P2 and P3. Indeed, the loss of introns compared to the ancestral structure of the gene and the trace of a polyadenylation tail are evidence to conclude that a pseudogene arose by retrotransposition.

Pseudogene analysis at the nucleotide level using phylogenetic tools

The mutation analysis and its interpretation in an evolutionary perspective allow the analysis of the pseudogenization process.

As a first step, all the mutations occurred in the pseudogenic sequence must be identified. This can be done by performing a comparative analysis between the nucleotide sequence of the pseudogene under study and that of the corresponding functional gene(s) (Fig. 6). The functional reference sequence may be selected among the available orthologs. For a better accuracy, it is also possible to reconstruct the ancestral sequence and use it as reference [START_REF] Dainat | GLADX: An automated approach to analyze the lineage-specific loss and pseudogenization of genes[END_REF][START_REF] Paten | Genome-wide nucleotide-level mammalian ancestor reconstruction[END_REF]. This is particularly suited when the divergence between the contemporary functional reference sequence and the pseudogene sequence is high.

Secondly, when all mutations have been documented, the use of a parsimony method over the defined species tree allows revealing their evolutionary history. We can now distinguish between the pseudogenization events that are shared and the independents ones. As exemplified in Figure 6, we can thus identify the first mutation that led to the pseudogenization and the history of mutations that accumulated ever since.

Automation

With the increasing interest in the field of pseudogene function and evolution, the phylogeny methods become essential. The automation of these methods is an important issue as well, but it is still in its infancy. Many difficulties exist, such as the construction of reliable phylogenies, the phylogenetic analysis, the collection of evidence from heterogeneous data, and the ability to merge all the information available to deduce the best explanation among the diversity of possibilities. To our knowledge, the only available automated method using phylogeny for pseudogene analysis is GLADX (Gene Loss Analyzer DAGOBAH eXtension) [START_REF] Dainat | GLADX: An automated approach to analyze the lineage-specific loss and pseudogenization of genes[END_REF]. It follows the outline described in this chapter. It relies on many analytical tools (e.g. Blast, phylogeny, protein prediction, reconstruction of ancestral sequences), and integrates heterogeneous data. Based on an expert system [START_REF] Gouret | Integration of evolutionary biology concepts for functional annotation and automation of complex research in evolution: The multi-agent software system DAGOBAH[END_REF], GLADX is in fact designed to be as close as possible to human expertise and, compatibly with the errors that can be present in the databases and with the limits of the tools used, it is able to provide reliable results and interpretations. However, GLADX has been designed to study specifically unitary pseudogenes and works exclusively with the Ensembl database. Furthermore, it has not been tested on recent databases and can be outdated. Finally, it is not maintained. The complexity of the framework in which it has been developed, as well as the uncommon programming language used (PROLOG), didn't help to find contributors to maintain it. Using new pipeline frameworks such as Snakemake [START_REF] Köster | Snakemake-a scalable bioinformatics workflow engine[END_REF] or Nextflow [START_REF] Di Tommaso | Nextflow enables reproducible computational workflows[END_REF] might simplify the automation of such analysis, probably going even beyond.

Closing remarks

The most accurate way to detect and analyse pseudogenes is by using phylogenetic approaches. The integration of evolutionary context provides unrivalled results for the understanding of pseudogenes and the underlying pseudogenization phenomenon. Unfortunately, the development of tools to automate the detection and analysis of pseudogenes using phylogenetic approaches is still required. Such tools should be able to analyse the three main types of pseudogenes: processed, unprocessed and unitary pseudogenes.

Notes 1.

When a phylogenetic tree needs to be built, the choice of the species and the time occurring between speciation events need to be chosen wisely, in order to have a data set that ensures high accuracy. An example is reported in Figure 7A. If the aim of the study is to attribute a date to the loss of a given gene observed in M. musculus, the use of H. sapiens and M. domestica as other species allows estimating that the loss occurred during a period of 100 million years. However, adding one more species such as R. norvegicus enables to deduce that the gene was lost during a more restricted period of 40 Million years (Fig. 7B). If the gene is absent in R. norvegicus as well, it can be deduced that it was lost after the Eutheria speciation and before the Rodent speciation during a period of 60 Million years. The loss might have occurred two times independently after the Rodent speciation (once in the R. norvegicus lineage and once in the M. musculus lineage), but this scenario would be less parsimonious and is excluded. 2. The choice of species must also consider the annotation quality of their genomes, because the more the genome is badly sequenced, the higher is the probability that a gene or pseudogene is present in a non-sequenced DNA fragment. Thus, plenty of missing genes and pseudogenes might in fact be sequencing artefacts and be actually present in the genome. The chance of finding the gene sequence of interest can be increased using EST databases and the absence or the pseudogenization of a gene can be double-checked using a resequencing process. In any case, a compromise has to be found between the choice of species and the quality of the annotation of their genomes. 3. When a pseudogene is ancient, its study should be performed at the amino acidic level, where the sequence can be more conserved due to synonymous mutations (Fig. 8). The same applies to a recent pseudogene, when the only known functional gene exists in a very divergent species. In all the other cases, it is possible to study the pseudogene at the nucleotide level. 4. To build the phylogeny, firstly a blast should be performed using a query sequence against different genomes in order to find similar genes. Using these sequences in FASTA format, a multiple sequence alignment should be then performed using tools such as ClustalW, Muscle or T-coffee [START_REF] Jungreis | Evolutionary dynamics of abundant stop codon readthrough[END_REF][START_REF] Stark | Discovery of functional elements in 12 Drosophila genomes using evolutionary signatures[END_REF][START_REF] Loughran | Evidence of efficient stop codon readthrough in four mammalian genes[END_REF]. To increase the accuracy of the phylogeny, the alignment can be cleaned in order to remove sites, sequences and regions that are spurious or poorly aligned [START_REF] Dainat | GLADX: An automated approach to analyze the lineage-specific loss and pseudogenization of genes[END_REF][START_REF] Paten | Genome-wide nucleotide-level mammalian ancestor reconstruction[END_REF][START_REF] Gouret | Integration of evolutionary biology concepts for functional annotation and automation of complex research in evolution: The multi-agent software system DAGOBAH[END_REF]. Once the alignment is ready, the phylogeny can be built using software such as Phyml, Figenix, RaxML [START_REF] Paganini | Reliable Phylogenetic Trees Building: A New Web Interface for FIGENIX[END_REF][START_REF] Stamatakis | RAxML-VI-HPC: maximum likelihood-based phylogenetic analyses with thousands of taxa and mixed models[END_REF][START_REF] Guindon | New algorithms and methods to estimate maximum-likelihood phylogenies: assessing the performance of PhyML 3.0[END_REF]. Note that Figenix integrates these steps and can automatically build a phylogeny from a query sequence.

5. An alternative method to detect gene losses is the use of clustering algorithms. Indeed, after the clustering of sequences of different genomes, the losses can be determined by a parsimony analysis of the genes contained in the clusters compared to the species tree of the species used. However, the use of the clustering method by itself can lead to artefactual groups. For large-scale studies, this approach has the advantage compared to a phylogenetic approach to be relatively fast. Thus, it can be used as first quick filter to be followed by a deepened analysis using phylogeny. The IODA database contains losses that have been detected in Chordate using this double approach [START_REF] Levasseur | The chordate proteome history database[END_REF]. 6. Two approaches are available to build this new phylogeny: i) to use the sequence alignment employed to build the previous phylogeny (see Note 4) after having added the putative pseudogenic sequence; ii) to build a new alignment with the sequences obtained by Blast using the putative pseudogene as a query. Using the first approach, if the putative pseudogene comes from a distant family, it might be placed as an external sequence. Thus, no information about its evolution or its function can be deduced. On the contrary, using the second approach, the new sequences collected for this phylogeny can provide new information. 7. Depending on the method used to build the new phylogeny, the gene initially used as query in the Blast step may be missing. Therefore, using only this gene to test orthology is risky. 8. If the sequence investigated is too degenerated to be kept in the process (too short or too divergent from the other present sequences), the orthology cannot be determined by phylogeny. 9. Often, the first Blast hits found are not pseudogenes, but paralogous or co-orthologous functional genes (since pseudogenes evolve more quickly than genes, their similarity is often lower than those of paralogous genes), while the hits that correspond to pseudogenes do not appear among the best. Moreover, the Blast search can detect several pseudogenes. Consequently, depending on what is sought for (only one pseudogene or as many as possible), it is necessary to study a bigger or smaller fraction of the Blast hits found. 10. As annotation errors or mis-annotations may exist in databases, the orthologous sequence retrieved may be a pseudogene annotated as functional gene or, on the contrary, a functional gene annotated as pseudogene. These cases are rare and difficult to assess therefore they require deeper investigation. The suspicious to face up such case can be raised when the pattern of the loss/pseudogenization events detected in the results are not parsimonious [START_REF] Dainat | GLADX: An automated approach to analyze the lineage-specific loss and pseudogenization of genes[END_REF]. 11. The checking may be done verifying the annotations present at the genome location where the sequence has been found. 12. Even if disruptive mutations are found, this does not mean that the identified pseudogene is non-functional. Pseudogenes have for a long time been considered without functional relevance, but numerous studies have indeed shown the transcriptional activity of some of them, and their involvement in biological processes [START_REF] Chan | Transcribed pseudogene ψPPM1K generates endogenous siRNA to suppress oncogenic cell growth in hepatocellular carcinoma[END_REF][START_REF] Hirotsune | An expressed pseudogene regulates the messenger-RNA stability of its homologous coding gene[END_REF][START_REF] Wen | Pseudogenes are not pseudo any more[END_REF][START_REF] Pink | Pseudogenes: pseudo-functional or key regulators in health and disease?[END_REF]. However, the transcriptional or functional activity may be tough to detect and requires further analyses. Furthermore, when a functional activity is revealed, it remains to be determined whether this activity is a residual activity of the parental gene or if the pseudogene can be assimilated to a new gene having a new function. 13. Depending on the chosen approach and on the subsidiary information desired, it may be preferable to perform the verification that the retrieved sequence is indeed a pseudogene before performing the orthology verification by phylogeny. It is a considerable time optimization to avoid building a phylogeny to study a sequence that finally would turn out to be an already annotated gene.

  It maintains the original function; this process is known as gene conservation.

o o It acquires a new function by neo-functionalization. o It shares the original function with the gene it originated from by subfunctionalization.

Acknowledgements

Thanks to Philippe Gouret and Julien Paganini for their help on the development of the strategy to detect gene losses and unitary pseudogenes.