Vincent Henry
email: vincent.henry@inria.fr-olivier.colliot@upmc.fr

Ivan Moszer

Olivier Dameron

Marie-Claude Potier

Martin Hofmann- Apitius

Olivier Colliot

Converting Alzheimer's disease map into a heavyweight ontology: a formal network to integrate data

Keywords: Alzheimer's disease, ontology, disease map, model consistency

Alzheimer's disease (AD) pathophysiology is still imperfectly understood and current paradigms have not led to curative outcome. Omics technologies offer great promises for improving our understanding and generating new hypotheses. However, integration and interpretation of such data pose major challenges, calling for adequate knowledge models. AlzPathway is a disease map that gives a detailed and broad account of AD pathophysiology. However, AlzPathway lacks formalism, which can lead to ambiguity and misinterpretation. Ontologies are an adequate framework to overcome this limitation, through their axiomatic definitions and logical reasoning properties. We introduce the AD Map Ontology (ADMO), an ontological upper model based on systems biology terms. We then propose to convert AlzPathway into an ontology and to integrate it into ADMO. We demonstrate that it allows one to deal with issues related to redundancy, naming, consistency, process classification and pathway relationships. Further, it opens opportunities to expand the model using elements from other resources, such as generic pathways from Reactome or clinical features contained in the ADO (AD Ontology). A version of ADMO is freely available at http://bioportal.bioontology.org/ontologies/ADMO.

Introduction

Alzheimer's disease (AD) is a progressive neurodegenerative disorder of the brain that was first described in 1906. The intense activity of AD research constantly generates new data and knowledge on AD-specific molecular and cellular processes (a Medline search for "Alzheimer disease" results in over 135,000 articles, as of June 30, 2018). However, the complexity of AD pathophysiology is still imperfectly understood [START_REF]2018 Alzheimer's disease facts and figures[END_REF]. These 110 years of efforts have essentially resulted in one dominant paradigm to underline the causes of AD: the amyloid cascade [START_REF] Golde | Anti-aβ therapeutics in Alzheimer's disease: the need for a paradigm shift[END_REF]. Therapeutics targeting this pathway failed to lead to curative outcome for humans, strongly suggesting the need for alternative hypotheses about AD etiology. Since the turn of the century, omics technologies lead to a more comprehensive characterization of biological systems and diseases. The production of omics data in AD research opens promising perspectives to identify alternatives to the amyloid cascade paradigm. The current challenge is thus to integrate these data in an appropriate way, in order to propose new hypotheses and models about AD pathophysiology.

Systems medicine disease maps (DM) provide curated and integrated knowledge on pathophysiology of disorders at the molecular and phenotypic levels, which is adapted to the diversity of omics measurements [START_REF] Ogishima | AlzPathway, an Updated Map of Curated Signaling Pathways: Towards Deciphering Alzheimer's Disease Pathogenesis[END_REF]- [START_REF] Fujita | Integrating Pathways of Parkinson's Disease in a Molecular Interaction Map[END_REF]. Based on a systemic approach, they describe all biological physical entities (i.e. gene, mRNA, protein, metabolite) in their different states (e.g. phosphorylated protein, molecular complex, degraded molecule) and the interactions between them [START_REF] Kitano | Using process diagrams for the graphical representation of biological networks[END_REF]. Their relations are represented as biochemical reactions organized in pathways, which encode the transition between participants' states as processes. AlzPathway is a DM developed for AD [START_REF] Ogishima | AlzPathway, an Updated Map of Curated Signaling Pathways: Towards Deciphering Alzheimer's Disease Pathogenesis[END_REF]. It describes 1,347 biological physical entities, 129 phenotypes, 1,070 biochemical reactions and 26 pathways. The information contained in DM is stored in syntactic formats developed for systems biology: the Systems Biology Graphical Notation (SBGN) [START_REF] Mi | Systems Biology Graphical Notation: Activity Flow language Level 1 Version 1.2[END_REF] and the Systems Biology Markup Language (SBML) [START_REF] Smith | SBML Level 3 package: Hierarchical Model Composition, Version 1 Release 3[END_REF]. While syntactic formats are able to index information, they are not expressive enough to define explicit relationships and formal descriptions, leading to possible ambiguities and misinterpretations. For AlzPathway, this defect in expressiveness results in the lack of formalism and thus of: a) hierarchy and disjunction between species (e.g. between "Protein" and "Truncated Protein" or between "Protein" and "RNA", respectively), b) formal definition of entities (such as phenotypes), c) formal relationships between reactions and pathways (that are missing or are managed as cell compartments), d) uniformity of entities' naming (e.g. complexes that are labelled by their molecular components or by a common name) and e) consistency between reactions and their participants (e.g. translation of genes instead of transcripts).

Compared to syntactic formats, the Web Ontology Language (OWL), a semantic format used in ontologies, has higher expressiveness [START_REF] Schaffert | A semantic wiki for collaborative knowledge formation[END_REF] and was designed to support integration. It is thus a good candidate to overcome the previous limitations.

An ontology is an explicit specification of a set of concepts and their relationships represented in a knowledge graph in semantic format. Ontologies provide a formal naming and definition of the types (i.e. the classes), properties, and interrelationships between entities that exist for a particular domain. Moreover, knowledge and data managed by an ontology benefit from its logical semantics and axiomatic properties (e.g. subsumption, disjunction, cardinality), which supports automatic control of consistency, automated enrichment of knowledge properties and complex query abilities [START_REF] Mizoguchi | Tutorial on ontological engineering[END_REF].

The Alzheimer's Disease Ontology (ADO) [START_REF] Malhotra | ADO: a disease ontology representing the domain knowledge specific to Alzheimer's disease[END_REF] is the first ontology specific to the AD domain. ADO organizes information describing clinical, experimental and molecular features in OWL format. However, the description of the biological systems of ADO is less specific than that of AlzPathway.

Considering that 1) semantic formats can embed syntactic information, 2) DM provide an integrative view adapted to omics data management and 3) an ontological model is appropriate to finely manage data, the conversion of AlzPathway into a formal ontology would bring several assets, including an efficient integration of biomedical data for AD research, interconnection with ADO and an increased satisfiability of the resources.

We propose the Alzheimer Disease Map Ontology (ADMO), an ontological upper model able to embed the AlzPathway DM. Section 2 is devoted to the description of the ADMO model. In Section 3, we describe a method to convert AlzPathway in OWL and how ADMO can manage the converted AlzPathway and automatically enhance its formalism. Section 4 presents elements of discussion and perspectives.

Ontological upper model: Alzheimer Disease Map Ontology

The initial definition of an ontological model aims to design a knowledge graph that will drive its content. In a formal ontology, the relationships are not only links between classes, but also constraints that are inherited by all their descendants (subclasses). Thus, the choices of axioms that support high level classes and their properties are key elements for the utility of the model. The Systems Biology Ontology (SBO) [START_REF] Courtot | Controlled vocabularies and semantics in systems biology[END_REF] is a terminology that provides a set of classes commonly used to index information in SBML format. These classes conceptualize biological entities at an adequate level of genericity and accuracy that supports a wide coverage with few classes and enough discrimination. We selected a set of 54 SBO terms from "process" or "material entity" for reactions and molecules as a first resource of subclasses of processes and participants, respectively. The modified Edinburg Pathway Notation (mEPN) [START_REF] Freeman | The mEPN scheme: an intuitive and flexible graphical system for rendering biological pathways[END_REF] is another syntactic format based on systems approach. Its components provide a refined set of molecular states that complete the SBO class set. Following class selection from SBO and mEPN, we designed a class hierarchy between them. We systematically added disjointness constraints between the generic sibling subclasses of participants in order to ensure that process participants belong to only one set (e.g. a gene cannot be a protein and reciprocally). We did not apply the same rule to the processes' subclasses as a reaction may refer to different processes (e.g. a transfer is an addition and a removal).

Properties consistent with a systems approach (i.e. part_of, component_of, compo-nent_process_of, has_participant, has_input, has_output, has_active_participant, de-rives_from and their respective inverse properties) were defined from the upper-level Relation Ontology (RO) [START_REF] Smith | Relations in biomedical ontologies[END_REF]. Then, we formally defined our set of classes with these properties and cardinalities to link processes and participants with description logic in SHIQ expressivity (e.g. a transcription has at least one gene as input and has at least one mRNA as output; a protein complex formation has at least two proteins as input and has at least one protein complex as output).

The design of the ADMO upper ontological model based on SBO, mEPN, RO and personal addition resulted in 140 classes (42 processes' subclasses and 83 participants subclasses) and 11 properties formally defined by 188 logical axioms in description logic (Fig. 1). This model is based on a simple pattern as our knowledge graph involves only three types of properties: 1) the is_a (subclass_of) standard property, 2) the has_part standard property and its sub-properties has_component and has_compo-nent_process and 3) the has_participant property and its sub-properties has_input, has_output and has_active_participant.

AlzPathway conversion and integration into ADMO

AlzPathway elements were extracted and stored in a structured table using home-made Python scripts. In this table, each biological entity was indexed by one of the high-level participants' subclasses of ADMO and all processes were in correspondence with their participants. The table also contains class annotations such as the AlzPathway identifier (ID), and IDs from other knowledge bases such as UniProt [START_REF]UniProt: the universal protein knowledgebase[END_REF] for participants and KEGG [START_REF] Kanehisa | KEGG: new perspectives on genomes, pathways, diseases and drugs[END_REF] for processes. The table is structured to integrate component information for multiplex entities (e.g. protein complex) and location information for the process (e.g. cell type or cell part). The table was then manually curated as described below.

In AlzPathway, native and modified proteins (e.g. phosphorylated or activated) may have the same label and same Id. In order to specify these different states, we added a suffix to modified protein labels (e.g. "_P" or "_a" for phosphorylated or activated, respectively).

In AlzPathway, phenotypes are participants. But several of them are named with a process name, pathway label or molecule type (e.g. microglial activation, apoptosis or cytokines, respectively). In order to deal with these ambiguities, 26 phenotypes were reclassified as molecules (e.g. cytokine) or cellular components (e.g. membrane) and 14 names that referred to processes or pathways were changed into processes' participant names (e.g. apoptosis became apoptotic signal). In addition, 5 phenotypes that were named with a relevant pathway name (e.g. apoptosis) were added to the initial set of the 26 AlzPathway's pathways.

AlzPathway only describes a subset of genes, mRNA and proteins. As omics technology can capture data at the genome, transcriptome or proteome levels, we added missing information in order to complete some correspondences between genes and gene products. This resulted in the addition of 406 genes, 415 mRNA and 194 proteins and protein complex states.

Then, using the ontology editor Protégé, the content of the structured table was imported into ADMO using the Protégé Cellfie plugin. Entities information were integrated as subclasses of ADMO participants classes. During the integration, we also added a new property has_template (sub-property of derives_from) to formally link a gene to its related mRNA and a mRNA to its related protein. Reactions were integrated as independent subclasses of the "process" class. Then, automated reasoning was used to classify them as subclasses of the ADMO upper model process classes depending on their formal definition (see Fig. 2a*). The 1,065 inferred subclass_of axioms corresponding to this refined classification of processes were then edited. During their import, process classes from AlzPathway were formally linked to their respective location through the RO property: occurs_in.

While AlzPathway does not formally link pathways and their related biochemical reactions, pathways were manually imported. For each pathway, a class "reaction involved in pathway x" was created and defined both as "reaction that has_participant the molecules of interest in x" and "component_process_of pathway x". For example, the class "reaction involved in WNT signaling pathway" has_participant "WNT" and is a component_process_of "WNT signaling pathway". Then, using automated reasoning, all reactions having participants involved in pathway x were classified as subclasses of "component_process_of pathway x" classes and were linked to the pathway with the component_process_of property by subsumption. For example, "SFRP-WNT association" is automatically classified as subclass of "reaction involved in WNT signaling pathway" (see Fig. 2b*) and inherits from its properties component_ process_of "WNT signaling pathway" (see Fig. 2b**). The 355 inferred subclass_of axioms corresponding to reactions involved in one of the 22 pathways were then edited. As a result, ADMO embeds AlzPathway in a consistent network containing 2,132 classes (2,175 disjoint participants, including 88 phenotypes or signals, 1,038 disjoint processes and 22 pathways) in relation with 10,964 logical axioms before and 12,373 logical axioms after automated reasoning, respectively. Specific efforts were dedicated to the design of classes hierarchy and formal definition with description logic axioms, leading to explicit relations between processes and biological entities. These axioms were inherited by classes imported from AlzPathway, resulting in the formal and precise description of the elements of AD pathophysiology. Thus, following automated reasoning, only 21 out of 643 AlzPathway's reactions generically classified as "transition" or "unknown transition" remained unaffected to a specific process of the ADMO upper model, such as metabolic reaction, phosphorylation or activation. Moreover, misaffected processes were consistently affected to a specific process (e.g. translation instead of transcription). In addition, 355 reactions were formally defined as subprocess of pathways of interest.

The conversion of AlzPathway also benefits from the ADMO simple pattern of relationships (Fig. 3) in which new properties were added: the derives_from property that links a modified protein to its native form, the has_template property that links a native gene product to its mRNA and gene, and the occurs_in property that links a process to its cellular location. Fig. 3. Alzheimer Disease Map Ontology (ADMO) pattern (A) and application to AlzPathway (B). AlzPathway classes (B; illustrated for the SFRP-WNT association process ant its participants) are now subclasses of ADMO classes (A). Each class of AlzPathway may be instantiated by the corresponding entities as individuals. Then, entities can be related to different objects in an RDF schema such as patients and experiments, or more specifically to values such as SNP for genes, relative expression for mRNA, and concentration for proteins.

Discussion

We proposed the ADMO ontological model in order to manage the conversion and integration of AlzPathway in OWL format. By converting AlzPathway into an OWL ontology, we increased its formalism. All entities are now formally defined and interconnected within a consistent network. While AlzPathway contained several ambiguities, our efforts on formalism at a semantic level for phenotypes and description logic in ADMO classes allowed us to solve inconsistencies. Moreover, the combination of SBO and mEPN provided a more precise specification of processes and biological entity states within the system compared to SBML or SBGN, which was beneficial for the specification of AlzPathway reactions following its import into ADMO. Unlike DM, ontologies are not adapted for graphical visualization but present a higher flexibility to integrate new elements in the knowledge graph, as we did by adding 865 genes and mRNA. Moreover, during the conversion step, AlzPathway's internal IDs were retained as class annotations, allowing interoperability between the initial and converted AlzPathway. Taking advantage of the knowledge graph and its semantic links, the ID information are retrievable from a derived molecule to its native form following the derives_from or has_component properties that link each of these classes.

Furthermore, the increased formalism requires to assert a participant as subclass of the most representative class and thus, clarifies the status of the entities. In several standard bioinformatics knowledge resources (e.g. UnitProt [START_REF]UniProt: the universal protein knowledgebase[END_REF], KEGG [START_REF] Kanehisa | KEGG: new perspectives on genomes, pathways, diseases and drugs[END_REF]), a same ID refers to a gene or a protein and in fine to a set of information, such as gene, interaction, regulation and post translation modification (PTM), which are thus not specifically discriminated. However, current omics technologies are able to generate data focused on specific elements of the systems (gene mutation, relative gene expression, protein concentration…). This is underexploited by standard resources. Based on DM approaches, we provided an ontology that a) represents the complexity of a system such as AD pathophysiology and b) is designed to specifically integrate each type of omics data as an instance of the explicit corresponding class.

The next possible step is to instantiate the model with biomedical omics data. To this end, the Resource Description Framework (RDF) semantic format is appropriate as it was specifically designed for representing a knowledge graph as a set of triples containing directed edges (semantic predicates). Different RDF schemas were already developed in the field of molecular biology (BioPax [START_REF] Demir | BioPAX -A community standard for pathway data sharing[END_REF]) or more specifically for AD biomedical research (neuroRDF [START_REF] Iyappan | NeuroRDF: semantic integration of highly curated data to prioritize biomarker candidates in Alzheimer's disease[END_REF]). The Global Data Sharing in Alzheimer Disease Research initiative [START_REF] Ashish | Global Data Sharing in Alzheimer Disease Research[END_REF] is also a relevant resource to help find appropriate predicates to enrich RDF schemas and refine subject information (age, gender, clinical visit…). Depending on the need of a given study, users may design RDF schemas with their own predicates of interest. Then, this RDF schema can be integrated in our ontology by adding data as instances of its corresponding specific classes (Fig. 3B). Therefore, instantiation opens perspectives for complex querying, both richer and more precise than indexing.

DM are based on systems biology approaches, allowing one to take each part of the system into consideration. Our ontology goes one step further by formally defining the different elements of the system and linking them with the biochemical reaction and pathway levels. Here, we relied on AlzPathway, but additional resources could be used, such as Reactome [START_REF] Fabregat | The Reactome Pathway Knowledgebase[END_REF] which provides a wide range of generic curated human biochemical reactions and pathways. Our ADMO upper ontological model provides an interesting framework to embed generic resources and thus harmonize AlzPathway and those resources. By converting and integrating AlzPathway in OWL format, the resulting ontology is ready to be connected with ADO and its clinical knowledge description. Owing to its specificity on biochemical reactions, an interoperable and formal version of AlzPathway should be a relevant complement to ADO. This offers new avenues for increasing the scale of representation of AD pathophysiology in our framework. In the same way, the genericity of processes and participants described in the ADMO upper model opens the perspective to harmonize specific DM from different neurodegenerative disorders such as the Parkinson's disease map [START_REF] Fujita | Integrating Pathways of Parkinson's Disease in a Molecular Interaction Map[END_REF] and others.

Fig. 1 .

 1 Fig. 1. Alzheimer Disease Map Ontology model design. Classes were extracted from the Systems Biology Ontology (SBO) and the modified Edinburg Pathway Notation (mEPN) into Protégé. Classes were hierarchized as subclasses of process (A) or participant (B). Using properties from the Relation Ontology (RO), classes were formally defined in description logic, as illustrated in the case of transcription (C) and protein complex formation (D) processes.

Fig. 2 .

 2 Fig.2. Example of automated reasoning on Protégé. Asserted axioms are in uncoloured lines and inferred axioms are highlighted in yellow. Following automated reasoning, SFRP-WNT heterodimer association is classified as subclass of "protein complex formation" (a*) and of "reaction involved in WNT signaling pathway" classes (b*), thus it inherits of the component_pro-cess_of "WNT_signaling pathway" property (b**).

Acknowledgements. The research leading to these results has received funding from the program "Investissements d'avenir" ANR-10-IAIHU-06 (Agence Nationale de la Recherche-10-IA Institut Hospitalo-Universitaire-6) and from the Inria Project Lab Program (project Neuromarkers).


Benjamin Rotenberg

Ignacio Pagonabarraga

Daan Frenkel

Coarse-grained simulations of charge, current and flow in heterogeneous media

Keywords:

We present a coarse-grained simulation method for complex charged systems. This mesoscopic model couples a hydrodynamic description to a free energy functional accounting for the interactions between solvent(s) and charged solutes. It is implemented in a hybrid lattice-based algorithm, whereby the evolution of the overall mass and momentum is taken care of via a Lattice Boltzmann scheme, whereas the composition and ionic concentrations are updated using the link-flux method. Several applications illustrate the power of this coarse-grained model for charged heterogeneous media : The transport of charged tracers in charged porous media, the deformation of an oil droplet in water under the effect of an applied electric field, and the distribution of ions at an oil-water interface as a function of their affinity for both solvents.

INTRODUCTION

The study of complex fluids and heterogeneous materials offers significant challenges because of the wide range of relevant length and time scales involved. Flow in porous media, the equilibrium and kinetic properties of membranes, oil/water mixtures or the electrokinetics of drops and colloids are examples of situations where the system evolves on scales much larger than molecular ones. Even the equilibrium structures characteristic of these systems are orders of magnitude larger than the molecular scale. Thus a simple atomistic description is not suitable and the use of soft coarse-grained potentials has become a standard alternative [START_REF] Karttunen | Novel Methods in Soft Matter Simulations[END_REF]. Under the action of external forces, the situation in these heterogeneous materials is even more involved, because of the hydrodynamic coupling between solutes and solvent. Hence, also the dynamic behaviour needs to be accounted for within a consistent coarse-grained approach.

A variety of strategies have been proposed to face these challenges. The combination of soft potentials and local thermostats which conserve momentum locally, as Dissipative Particle Dynamics [2] or the Lowe-Andersen thermostat [3], has allowed to reach hydrodynamic scales while keeping track of some microscopic details. Ideas from kinetic theory have also lead to flexible tools. Stochastic Rotation Dynamics couples molecular solutes to a coarse-grained solvent which recovers hydrodynamic behaviour [4]. Lattice Boltzmann (LB) [START_REF] Succi | The Lattice Boltzmann Equation for Fluid Dynamics and Beyond[END_REF], a method evolved from lat-tice gases to describe the hydrodynamics of fluids at long scales, has been extended to account for complex fluids. In particular, the combination of free energy based models coupled to a hydrodynamic description was introduced by Yeomans et al. [START_REF] Swift | [END_REF][7][8] for simulations of non-ideal fluids and binary mixtures and later extended for other systems such as binary mixtures with surfactants [9][10][11] liquid crystals [12][13][14][15], ternary mixtures [16] or active fluids [17]. LB methods are particularly well suited for hydrodynamic simulations of fluids, especially in complex media, for it is in principle easy to parallelize the codes and to implement boundary conditions at solid/fluid interfaces [18]. This allowed the simulation of binary fluids in porous media [19], colloids at an oil-water interface [20][21][22] and suspensions of charged colloids [23].

Standard LB schemes for complex fluids often lead to spurious fluxes across solid-fluid boundaries or at liquid-liquid interfaces [24]. This can become catastrophic when one deals with charged solutes, as such fluxes could result in a progressive breakdown of electroneutrality. To overcome this difficulty, Capuani et al. [25,26] introduced the link-flux method to reproduce the convective-diffusive dynamics of charged species in an electrolyte. Within this hybrid scheme, the overall mass and momentum of the fluid are evolved using a LB algorithm, whereas the ionic concentrations are updated using the link-flux method (for diffusion and migration) in combination with an advection scheme (for convection). This algorithm satisfies detailed balance at steady-state and allows to rigorously cancel fluxes into the solid, even in pres-ence of moving boundaries.

We want to put forward a general scheme to deal with complex charged fluids making use of this hybrid strategy, by combining the LB treatment of the solvent with a continuum treatment of binary charged fluids. We show how starting from a general free energy functional it is possible to obtain a dynamic scheme consistent with the prescribed thermodynamics and how it can be coupled to the LB method for the fluid flow. Several applications illustrate the power of this coarse-grained model for charged heterogeneous media : The transport of charged tracers in charged porous media, the deformation of an oil droplet in water under the effect of an applied electric field and the distribution of ions at an oil-water interface as a function of their affinity for both solvents. We compare with analytic predictions in the appropriate limits to illustrate that it is possible to achieve quantitative control of the performance of the model.

FREE-ENERGY BASED MODEL OF NON-IDEAL MIXTURES

In this section, we present the description of a complex mixture of solvents and charged solutes in terms of mesoscopic variables, the associated thermodynamic properties and the hydrodynamic equations governing their dynamics.

Mesoscopic description

Our aim is to describe charged solutes in heterogeneous media involving solid-liquid or liquid-liquid (e.g. oil-water) interfaces. For not too concentrated solutions, we can treat the fluid as a continuum, whose state is characterized by its local mass density ρ(r, t) and concentrations in each solute ρ k (r, t). In the following we will always consider 1:1 electrolytes and denote the corresponding solutes by + and -for cations and anions, respectively. For oil-water mixtures, we also introduce the local composition

φ(r) = ρ o (r) -ρ w (r) ρ o (r) + ρ w (r) ∈ [-1; 1] (1)
where the w and o subscripts refer to water and oil, respectively. Immiscible fluids are characterized by regions where φ is almost constant (φ ∼ -1 in water and φ ∼ +1 in oil) separated by a "sharp" interface. The underlying assumption in the above description is to consider the fluid can be seen as locally homogeneous, although density, composition and solute concentrations can vary on a larger scale.

Thermodynamics

The thermodynamics of the system is determined by its free energy, expressed as the functional :

F [φ, ρ + , ρ -] = dr F V [φ(r), ρ + (r), ρ -(r)] (2)
where F V is a free energy density. Following Onuki [27][28][29] we separate the contributions F mix [φ] describing the immiscible solvents and F ions [φ, ρ + , ρ -] describing ions in a solvent of composition φ. The mixing contribution is chosen of the Landau-Ginzburg form :

F mix = dr - 1 2 Bφ 2 + 1 4 Bφ 4 + 1 2 K(∇φ) 2 . ( 3
)
The first two terms correspond to the bulk phase behaviour, with minima for φ = ±1, while the last reflects the cost of sustaining interfaces. The dimension of B is energy×length -3 and that of K is energy×length -1 . This standard choice gives at equilibrium a planar interface of the form φ(x) = tanh(x/ξ) with a width ξ = 2K/B and a surface tension σ = 8KB/9. The ionic contribution of the free energy consists of an ideal, a solvation and an electrostatic term :

F ions = dr α=± ρ α (r) [k B T (ln ρ α (r) -1) -µ α +V solv α (r) + z α e 2 ψ(r) (4)
where z ± = ±1 is the valency of the ions, µ α is a reference chemical potential and the electrostatic potential ψ is solution of the Poisson equation :

∇ • [ (r)∇ψ(r)] = -[ρ + (r) -ρ -(r)] e . (5)
The dielectric constant (r) depends on the local composition of the fluid. Although it could be a priori an intricate function of φ, it is reasonable to assume a linear relation (r

) = ¯ [1 -γφ(r)] with ¯ = w + o 2 the average dielectric constant and γ = w -o w + o ∈ [0; 1] the dielectric contrast.
The ionic solvation potential V solv ± accounts for the different solvation free energy in the two solvents. It is therefore natural to parametrize it as a function of the composition as V solv

± (r) = ∆µ ± 1+φ(r) 2
, where we have introduced for each ion the solvation free energy difference between water and oil ∆µ ± = µ o ± -µ w ± , also referred to as extraction or Gibbs transfer free energy in the electrochemistry community. As the free energy (2) does not exclude in principle values of the composition parameter outside the φ ∈ [-1, 1] range, we have considered that the physical properties ( , V solv ± ) for φ < -1 are that of water ( w , 0) and for φ > 1 that of oil ( o , ∆µ ± ).

Let us now briefly analyze the properties of this free energy in terms of the chemical potentials associated to ρ ± and φ and the corresponding thermodynamic forces acting on the fluid. The ionic chemical potentials are of the usual form :

µ ± = δF δρ ± = k B T ln ρ ± + V solv ± + z ± eψ , (6)
while the chemical potential corresponding to the solvent mixture reads :

µ φ = δF δφ = µ mix φ + µ solv φ + µ el φ . (7)
The first term is simply :

µ mix φ = -Bφ(r) + Bφ(r) 3 -K∇ 2 φ(r) . (8)
When the solvation free energy is taken as

V solv ± (r) = ∆µ ± 1+φ(r) 2
the second term is :

µ solv φ = [ρ + (r)∆µ + + ρ -(r)∆µ -] 2 . (9)
Finally, the electrostatic contribution is :

µ el φ = E(r) 2 2 δ δφ = - γ¯ 2 E(r) 2 . ( 10
)
Excess chemical potential gradients give rise to a thermodynamic force (per unit volume) that can be expressed as a pressure gradient from the Gibbs-Duhem equality :

f th V = -∇P = φ∇µ φ + ρ + ∇µ ex + + ρ -∇µ ex - ( 11
)
where the ex superscript refers to the excess chemical potentials. Each part of the chemical potentials (mixing, solvation and electrostatic) contribute to this force. For example, gradients of µ ex ± give rise to the force ρ el E, with E = -∇ψ the electric field, and gradients of µ el φ are the source of the dielectrophoretic force ( (r) -¯ )∇ E 2 2 . The latter drives oil-rich fluid elements (φ > 0, < ¯ ) towards region where E 2 is small and water-rich fluid elements towards region of higher E 2 (note that ∇µ el φ also generates a composition flux in addition to this force) and is particularly important in phenomena such as electrowetting. These electrostatic contributions to chemical potential gradients illustrate some new features captured by our free energy model compared to previous ones used in Lattice Boltzmann simulations of binary mixtures.

Hydrodynamics

The thermodynamic description of the system needs to be supplemented by a set of prescriptions for the dynamics. Overall mass conservation of the fluid implies :

∂ t ρ + ∇ • (ρu) = 0 ( 12
)
with u the local barycentric velocity of the fluid. Momentum conservation of the fluid and viscous dissipation then enter in the Navier-Stokes equation, which reads for and incompressible fluid :

∂ t u + u • ∇u = ν∇ 2 u + f V ρ ( 13
)
where ν = η/ρ is the kinematic viscosity and f V = f ext V +f th V is the sum of the externally applied force and thermodynamic force (11). Although we have not considered this in the following, it is in principle possible to introduce a composition-dependent viscosity [30].

The composition φ and ionic concentrations ρ ± also satisfy conservation laws :

∂ t φ + ∇ • (φu) = -∇ • j φ (14a) ∂ t ρ ± + ∇ • (ρ ± u) = -∇ • j ± (14b)
where we have introduced fluxes in the barycentric frame j φ and j ± . The latter are given by phenomenological equations, namely :

j ± = -D ± ρ ± ∇(βµ ± ) (15)
with D ± the ionic diffusion coefficients, and the Cahn-Hilliard equation for the composition :

j φ = -M ∇µ φ (16)
with M a mobility. The units of M differ from those of D ± and a diffusivity (in m 2 s -1 ) of the interface can be defined as D φ = M B(-1 + 3φ 2 ) with B from (3). Note that near the interface φ ∼ 0 and D φ < 0 : this "negative diffusion" maintains the composition jump at the interface.

Discussion

The above description of mixtures of solvent and ions and the particular choice of free energy functional are very similar to the ones adopted by Onuki [27][28][29]. The functional differs only on two points. First, we follow previous LB studies of binary mixture [18,20,31,32] and use the Landau-Ginzburg (LG) functional (3) instead of the Bragg-Williams one for numerical convenience. The combination of

LG and electrostatic free energies (without solvation terms) has also been used to investigate the wetting of a solid substrate by ionic solutions [33]. The second difference with Onuki's approach consists in using ρ el ψ for the electrostatic energy instead of (∇ψ) 2 . This more natural choice, consistent e.g. with the DFT work of van Roij et al. [34,35], doesn't require any approximations for the treatment of "image charges".

A major difference with both Onuki's and van Roij's work is that the present model captures not only the equilibrium states of the system, but also its dynamics. It shares many features with Dynamical Density Functional Theory [36][37][38][39], since it relies on an expression of fluxes proportional to gradients of chemical potentials. The free energy functional described above is relatively simple, for it neglects e.g. the effect of the finite size of the ions and correlations beyond the mean-field level. This description is perfectly valid for dilute solutions of 1:1 electrolytes and can be improved if more concentrated solutions or multivalent ions are considered. Moreover, it captures the presence of immiscible solvents and the (possibly asymmetric) affinity of the ions for one of them. The free energy model put forward in this paper can be seen as a simple limiting case of more elaborate free energy functionals. The essential difference with DDFT is that the hydrodynamic behaviour of the fluid is properly described (it satisfies the Navier-Stokes equation resulting from momentum conservation and viscous dissipation), whereas most DDFT studies consider a fluid at rest or mediating hydrodynamic interactions between large solutes via effective interactions (Oseen or Rotne-Prager tensors) [40]. This latter approach can be efficient for suspensions of solid particles, but is not valid a priori for liquid droplets in another liquid.

LATTICE SIMULATIONS

The coarse-grained model introduced in the previous section couples a hydrodynamic description of the fluid to a free energy based representation of its thermodynamic behaviour. The purpose of the present section is to introduce the computational methods used to solve numerically the coupled evolution equations for the composition, ionic concentrations and fluid velocity. The general strategy relies on the use of different lattice models, which are described here successively.

Hybrid lattice scheme

For non-ideal multicomponent fluids, the standard LB approach treats all species on the same footing and populations (see below) are associated to each component. This can become computationally expensive for more than two components. Nevertheless, this approach has recently been applied for a mixture of two solvents and two reactive solutes by Furtado et al. [41] to study the convective drop motion driven by non-linear kinetics. An alternative is to use a hybrid LB / finite elements approach, whereby the LB fluid is described by populations evolving as before, but the order parameters (e.g. the composition φ) are described by scalar fields evolving according to finite elements schemes [42,43]. As mentioned in the introduction, such methods may suffer from spurious fluxes which motivated the development of the linkflux method by Capuani et al. [25,26]. Here we generalize the link-flux approach to the more complex case of ions in a mixture of solvents. The overall mass and momentum are taken care of via a Lattice Boltzmann scheme, whereas the composition φ and ionic concentrations ρ ± are updated using the link-flux method. We now develop these two steps.

Lattice Boltzmann

Lattice Boltzmann (LB) is a well established method for hydrodynamic simulations based on kinetic theory [START_REF] Succi | The Lattice Boltzmann Equation for Fluid Dynamics and Beyond[END_REF]44]. The Boltzmann equation is a mesoscopic kinetic equation which determines the evolution of the probability density function f (r, v, t) of finding a fluid particle with a velocity v at position r and time t. The hydrodynamic fluid variables are derived as moments in velocity space of the distribution function starting with the fluid local density ρ = f dv and mass flux ρu = f vdv. Although the relaxation of the distribution function toward equilibrium is determined by a nonlinear collision operator, the Bhatnagar-Gross-Krook (BGK) model [45] shows that proper hydrodynamics can be already recovered from a linearized collision operator if it is isotropic and conserves mass and momentum.

LB tracks the dynamics of fluid populations, f i = f (r, c i , t), which evolve on the nodes r of a lattice of spacing ∆x moving to neighbouring nodes at finite time steps through a finite set of allowed velocities {c i } i∈ [START_REF] Karttunen | Novel Methods in Soft Matter Simulations[END_REF]Nmax] . The hydrodynamic variables are obtained as appropriate quadratures ρ = i w i f i and ρu = i w i f i c i , where the weights w i are associated to the chosen set of velocities. The particle distributions f i also relax at each time step toward a prescribed equilibrium distribution through a linear collision operator which conserves mass and momentum and ensures that the the solvent mass density ρ and velocity u follow the Navier-Stokes equations, ( 12) and ( 13), on distances larger than ∆x. The natural units in LB simulations are the lattice spacing ∆x and the time step ∆t. They are fixed by the properties of the system: In the following, the lattice spacing is chosen as a fraction of the Bjerrum length l B = e 2 /4π k B T (approximately 0.7 nm in water at room temperature), while the time step is determined by the value of the solvent viscosity.

The force per unit volume f V acting on the fluid also enters in the collision rule. The issue of computing the thermodynamic force (11) from the state (φ, ρ ± ) of the system will be addressed in the next section. For numerical accuracy and stability reasons, the fluid velocity must remain small compared to the sound velocity c s . This small Mach number limit implies that the forces are also small, i.e. that χ T f V ∆x 1, with χ T the fluid compressibility (χ T = 1/ρc 2 s for the LB fluid). Each term in the excess free energy density contributes to the force (11) and we can derive conditions accordingly. This leads for the F mix contribution (3) to χ T B ξ/∆x and similar requirements are obtained for the solvation and electrostatic ones.

Link-flux method

The composition of the fluid is characterized by the order parameters φ and ρ ± and evolves according to Eqs. ( 14), (15) and (16). The numerical solution of these equations is achived by the link-flux method. This method was introduced by Capuani et al. [25,26] in order to prevent spurious solute fluxes across solidfluid interfaces. It focusses on solute fluxes between lattice nodes rather than the amount of solute at each nodes. Integrating the conservation laws ( 14) over a volume corresponding to one lattice node and using Green's formula, we associate the variation of φ and ρ ± to the fluxes of j φ and j ± across the surface A 0 of the cell (for more details, see Ref. [25]). The latter can be separated into the contributions of each link between the considered node and all its neighbours :

ρ ± (r, t + ∆t) -ρ ± (r, t) ∆t ∆x 3 = -A 0 i j i ± (r) ( 17
)
where i refers to the discrete velocities and j i ± to the contribution of link i between r and r + c i ∆t to the outward flux of j ± . A similar expression can be written for the composition with link-fluxes j i φ . In order to ensure that the ions follow a Boltzmann distribution at equilibrium, we rewrite (15) as

j ± = -D ± e -βµ ex ± ∇[ρ ± e βµ ex ±
] and express the linkfluxes in the symmetrized form :

j i ± (r) = -d ± e -βµ ex ± (r) + e -βµ ex ± (r+c i ∆t) 2 × ρ ± (r + c i ∆t)e βµ ex ± (r+c i ∆t) -ρ ± (r)e βµ ex ± (r) ∆ i (18)
with d ± = D ± /A 0 and ∆ i = ||c i ∆t||. For links crossing solid-fluid boundaries (i.e. such that r + c i ∆t is a solid node) we enforce j i ± = 0 so that such links do not carry any solute into the solid. For the solvent order parameter φ we use the simpler form :

j i φ (r) = -m φ µ φ (r + c i ∆t) -µ φ (r) ∆ i (19)
with m φ = M/A 0 to recover the Cahn-Hilliard expression ( 16).

The link-flux algorithm just described takes care of the diffusive fluxes j φ and j ± . The advective fluxes φu and ρ ± u are treated in a separate step described in detail in ref. [25]: It consists in transferring particles according to the overlap between the considered cell (around a node) translated by u∆t and each of its neighbouring cells. Finally, the update of the composition also modifies the local force (11) acting on the fluid. The thermodynamic force acting on node r is the average of the forces on each link f V (r) = i w i F i (r)c i , with :

F i (r) = -φ(r) j i φ (r) m φ -k B T α=± j i α (r) d α - ρ α (r + c i ∆t) -ρ α (r) ∆ i (20)
For the ionic part we substract the ideal term ρ ± ∇µ± id = ∇ρ ± from the link flux, because only the excess (solvation and electrostatic) term ρ ± ∇µ ex ± contributes to the force. The computation of the chemical potentials µ ± and µ φ require the evaluation at the lattice nodes of the electrostatic potential ψ and its gradient E = -∇ψ, as well as ∇ 2 φ (see Eqs. ( 8) and ( 10)). The electrostatic potential is determined from the local charge and dielectric constant by solving the Poisson equation [START_REF] Succi | The Lattice Boltzmann Equation for Fluid Dynamics and Beyond[END_REF] with the Successive Over Relaxation (SOR) method [START_REF] Press | Numerical Recipes in C: The Art of Scientific Computing[END_REF], which also requires to compute gradients and Laplacian of φ and ψ. We have implemented a modified SOR algorithm with the additional term arising from the spatial variations of the permittivity. The following stencils are used to compute these differential operators :

∇ψ(r) = 1 ∆t i w i c 2 s [ ψ(r + c i ∆t) -ψ(r) ] c i (21a) ∇ 2 ψ(r) = 2 ∆t 2 i w i c 2 s [ ψ(r + c i ∆t) -ψ(r) ] (21b)
This choice was motivated by numerical accuracy considerations, but also by the fact that it is fully consistent with the discretization used for the LB part of the hybrid scheme.

CHARGED TRACERS IN CHARGED POROUS MEDIA

We now show a first application of the method introduced in the previous section in which no oil is present. To assess the behaviour of a forced, charged fluid, we consider a slit of width L with charged walls (charge density σ el < 0) and compensating counterions in the fluid as a simple representation of a porous medium. In equilibrium, the counterion con- where the amplitude satisfies u ref ≡ eE 0 /2πηl B . In Fig. 1 we compare the theoretical predictions and the computed values for the flow field for increasing charge of the solid walls. One can see that the method describes quantitatively the osmotic flow deep into the non-linear regime of electrostatic coupling (αL → π).

centration is ρ + (x) = α 2
Previous work has shown that other dynamic quantities, such as the dispersion of charged tracers by this flow, can be easily recovered [START_REF] Rotenberg | [END_REF].

A very useful quantity to analyze the diffusive dynamics in porous media is the time-dependent diffusion coefficient D(t), that can be measured by NMR [48]. At short times D(t) coincides with the molecular diffusion coefficient D. For neutral tracers, the short-time behaviour of D(t) reflects the geometry of the pores [49]:

D(t) D ∼ 1 -4 9 √ π S Vp
√ Dt with V p the pore volume, S the surface of the solid. The long time limit of D(t) is the effective diffusion coefficient :

D e = lim t→∞ D(t) (23)
which reflects the connectivity between pores throughout the medium. The ratio D e /D is often referred to as the inverse of the tortuosity. The effective tracer diffusion coefficient through charged porous materials is known to depend on its charge. For example, it has been observed in clays (negatively charged minerals) that the ratio D e /D for cations is larger than for neutral tracers, whereas that for anions is smaller [50][51][52].

This can be at least partly explained by the so-called Donnan effect : The concentration of co-ions (resp. counterions) in the pores of the material is smaller (resp. larger) than the concentration of salt in the reservoirs used to impose a concentration gradient to the sample, and this modifies the effective concentration gradient inside the sample. Even if we correct for this effect, we expect that ions of different charge will follow different pathways through the pores and this might influence the observed value of D e . In the framework of the proposed model, we can analyze D(t) as the integral of the tracer's velocity autocorrelation function (VACF) using the Moment Propagation [START_REF] Rotenberg | [END_REF][53][54][55]. It is worth emphasizing that the present analysis is possible because we are able to evaluate D(t) numerically. Although this has been done for neutral tracers by averaging over trajectories of explicit tracer particles [56], we found no such analysis for charged tracers.

System

We will analyze the diffusion of charged tracers in a porous medium consisting of a compact FCC lattice of charged spheres of radius R with a surface charge density σ el < 0 whose pores are saturated with an electrolyte solution of concentration ρ b . The void frac- corresponds to the exponential decay of the potential near a planar interface, for not too high a surface potential ψ S (compared to k B T /e and using the potential of the solution "far" from the surface as a reference). The latter depends on the surface charge density and the salt concentration in the medium.

Simulations were performed on a a 3 = (100∆x) 3 lattice, with spheres of radius R = a/2 √ 2 ∼ 35.4∆x. The lattice spacing is 2.5l B ∼ 1.75 nm, so that the O h (resp. T d ) pore size is ∼25.8 nm (resp. 13.9 nm). The charge density of the solid is σ el ∼ -0.04 e/nm 2 , and we considered salt concentrations ρ b corresponding to (κ b R) -1 ∈ [0.02; 0.57]. The molecular diffusion coefficients are D ± = 5 10 -2 (∆x 2 /∆t). The system is initialized with cations and anions distributed homogeneously and evolved until the equilibrium distribution is reached.

Results and discussion

The time-dependent diffusion coefficient D(t) for tracers of valency ±1 and 0 with the same molecular diffusion coefficient D are computed using the moment propagation method. The results show that the charge of the ion influences both the value of D e and the transient regime to reach this asymptotic value. D e /D is larger (resp. smaller) for cations (resp. anions) than for neutral tracers. This is in agreement with experimental observations (which also reflect the Donnan effect). D(t) also tends towards D e faster (resp. more slowly) for cations (resp. anions) than neutral tracers. This can be quantified by a characteristic time :

τ = ∞ 0 D(t) -D D e -D dt (24)
Fig. 2 shows the variation of tracer's diffusion and the relaxation time τ with salt concentration ρ b . The extension of the diffuse layer (approximately κ -1 b ) decreases with increasing ρ b . Both D e and τ for the charged tracers tend towards the values for neutral tracers at high ρ b . This effect has been observed in diffusion experiments (also reflecting the smaller Donnan effect in that case). To our knowledge, there is no experimental measurements of time-dependent diffusion coefficient for ions in clays. The variations of D e /D and τ with the tracer charge and salt concentrations reflect how tracers go from one cavity to another to explore the whole porosity. Smaller T d pores act as bottlenecks through which tracers must pass to go from one O h pore to another. In addition to this purely geometric (entropic) effect, the electrostatic potential distribution in the porosity also affects the motion of charged tracers. The electrostatic potential ψ(r) in the pore is always larger than the surface potential ψ S and increases with increasing distance from the surfaces. Thus ψ is larger near the center of O h pores (see the disconnected isopotential surfaces for a large ψ, in red in Fig. 2) than in T d pores and anions feel a repulsive electrostatic force when approaching the latter. This decreases the probability to go from one O h cavity to the next and consequently D e is smaller than for netral tracers. This also implies that it takes longer for anions to explore the volume accessible to them (although this volume is smaller than for neutral species that are not repelled from the surfaces) and the corresponding τ is larger.

As opposed to anions, cations accumulate near the surface. As can be seen in Fig. 3, the diffuse layer forms a continuous volume throughout the porosity and cations can follow preferential pathways along the surfaces. This surface diffusion mechanism is more FIG. 3: Isopotential curves for a FCC lattice of charged spheres where the porosity contains solvent, counterions and salt. The grey surface is an electric isopotential curve for a value close to the surface potential while the red surfaces correspond to a larger value of the electrostatic potential. Each lattice node in the solid phase is represented as a blue dot, while fluid nodes are not indicated for clarity.

efficient than for neutral cations, since exploring a smaller volume takes less time than exploring the whole O h pores. Hence the larger D e and smaller τ for cations. In both cases, there is a clear interplay between the geometric and electrostatic effects. At higher concentrations, the potential variations are screened and the effect of the above-mentioned mechanisms are less pronounced.

DIELECTRIC DROPLETS UNDER AN ELECTRIC FIELD

In this section, we apply our mesoscopic model to the deformation of an oil droplet in water under an applied electric field, in the absence of ions. Because of surface tension, the equilibrium shape of an oil droplet in water corresponds to the minimal interface area (a disk in 2D, a sphere in 3D) and there exists an excess (Laplace) pressure inside the drop : The final shape is governed by the balance between electrostatic and surface tension forces. For small applied fields the equilibrium shape is an ellipse in 2D (an ellipsoid in 3D) and the deformation is defined as D = (b -a)/(b + a) with b (resp. a) its large (resp. small) axis. In the small E limit, an analytical result for the deformation in the 2D case can be obtained following the lines of [57] for spherical droplets, with the result (see below):

P L = σ R d in 2D and 2σ R d in 3D,
D theor = 1 4 γ 2 (1 + γ) ¯ E 2 R d σ ( 25
)
where ¯ = ( w + o )/2.

System

We performed simulations of a two-dimensional oil droplet in water. Because of periodic boundary conditions in all directions, this corresponds to an array of infinite cylinders. The box size is N × N × 1 lattice points, with N = 50 or 100. The parameters entering in the free energy (3) are βB∆x 3 = 10 -3 and βK∆x = 3 10 -3 , giving a theoretical interface width of ξ = 2K/B ∼ 2.45∆x. and a surface tension βσ∆x 2 = 8KB/9 ∼ 1.63 10 -3 . This choice of parameters ensures that the interface is thin while remaining well resolved on the lattice. In particular, we checked that in the case of a planar interface (1D geometry) the simulated systems reproduces accurately the φ(x) = tanh(x/ξ) profile with the expected width. These parameters also fulfill the condition χ T B ξ/∆x. The interface mobility M is such that M B = 10 -2 ∆x 2 /∆t. After equilibrating the system, we turn on the electric field E by solving the Poisson equation under the condition of a potential drop -EL between both sides of the simulation box (to be consistent with the periodic boundary conditions). The electrostatic potential inside and outside the droplet is reported in Fig. 4 together with the correspongindg electric field lines, for a large dielectric contrast γ = 0.9 and an applied electric field βeE∆x = 10 -3 . The average dielectric constant ¯ = ( w + o )/2 is such that lB = βe 2 /4π¯ = 0.4∆x. Simulation results for two system sizes are shown : The smaller (b) corresponds to a box size L = 50∆x and a droplet radius R d = 14.1∆x, the larger (d) to L = 100∆x, R d = 28.2∆x. These results are compared to the analytical solution for an isolated cylinder : The simulation results for a system with periodic boundary conditions are in good agreement with the analytical solution for an isolated droplet. In particular, the field lines tend to bypass the less dielectric droplet and are thus closer to each other in the water phase. The field strength is smaller in the water phase, where the dielectric constant is higher. The analytical solution (26) corresponds to an infinitely thin interface (ξ = 0) with a discontinuous dielectric constant, whereas in simulations the latter varies smoothly over a distance ∼ ξ. The effect of the finite width is more pronounced for the smaller system (ξ/R d ∼ 0.17) and the agreement with ( 26) is excellent for the larger (ξ/R d ∼ 0.09). Although smaller ξ/R d can be reached by decreasing ξ, the latter should remain at least of a few ∆x for the variations of the composition φ to be well resolved on the lattice.

E(r) = (1 + γ)E 0 , r < R d E(r) = E 0 + γR 2 d ∇ E 0 • r r 2 , r > R d (26)

Results and discussion

Under the applied electric field, gradients of µ φ lead to an evolution of the composition φ, corresponding to a deformation of the droplet. The precise location of the interface is defined by the curve φ(r) = 0, which is fitted to an ellipse to obtain the deformation D. Results are reported in Fig. 5 as a function of the dielectric contrast γ for two droplet radii and two electric field strengths. For the simulation parameters used here, the largest deformation is < 10 -2 so that we are always in the small deformation limit assumed in Eq. ( 25). Simulation results are in excellent agreement with the analytical ones. This demonstrate that the simulation method faithfully reproduces the prediction of the continuous free energy model. The deviations observable only for large values of γ depend mainly on the radius R d and less significantly on the applied electric field E. This can be understood by analyzing the effect of the finite width on the energy balance leading to the final shape of the droplet and the theoretical deformation (25). In the case of a thin interface (ξ = 0) and assuming that the field at steady-state is still given by ( 26), we find that the electrostatic energy stored in the droplet when the field is turned on is (per unit length of the cylinder) :

drop 1 2 ( o E 2 in -w E 2 0 )dS = -1 2 ¯ (1 + γ)γ 2 E 2 0 S, with S = πR 2
d the section area of the droplet. This term tends to increase S by deforming the droplet and is balanced by the increase in surface energy σ2πR d D (per unit length of the cylinder). The corresponding deformation D is given by Eq. (25).

The main effect of the finite width ξ is that the smooth variation of at the interface leads to a O(ξR d ¯ (1 + γ)γ 2 E 2 0 ) correction to the electrostatic term, resulting in a smaller driving force for the deformation. This explains the field-independent O(ξ/R d ) deviation of the simulation results from D theor /(¯ E 2 R d /σ) observed in Fig. 5. Other possible field-dependent corrections do not seem to be observed.

The results shown in this section demonstrate that the lattice simulation scheme is able to reproduce quantitatively the behaviour of a system consisting of two immiscible dielectric fluids under an electric field. We have not considered here the possibility of altering the phase behaviour of the solvent mixture in an inhomogeneous electric field, but this phenomenon is also captured by the free energy model used for our simulations [58].

IONS AT AN OIL-WATER INTERFACE

In this last section, we further exploit the power of our coarse-grained description by considering ions at a planar oil-water interface. The two solvents may have different dielectric constants and the ions different affinities for the two solvents (∆µ ± = µ o ± -µ w ± = 0). The salt is partioned between the two phases and the salt concentration ratio is ρ o /ρ w = exp(-β∆µ av ) with ∆µ av = (∆µ + + ∆µ -)/2. In the case of asymmetric solvation (∆µ + = ∆µ -) there exists at equilibrium an electrostatic potential difference (Donnan potential) across the interface ψ D = ψ o -ψ w = (∆µ --∆µ + )/2e. In that case there is a charge separation over distances characterized by two different Debye lengths

κ -1 o,w = (2βe 2 ρ o,w / o,w ) -1/2
. For large values of the dielectric contrast or of ∆µ av , they can differ by several orders of magnitude.

Experimental investigations of such interfaces at the microscopic scale have revealed that the interplay between solvation and electrostatic forces at the interface is rather complex [59]. This is particularly true if there is a dielectric constrast between the two phases (γ = 0) resulting in image charges interactions. While Onuki accounted for them by an effective "image charge potential" in his free energy description of mixtures of oil, water and ions [27,28], van Roij et al. [35] suggested a Poisson-Boltzmann (PB) treatment for a thin planar interface introducing a shift s between the true solvent interface (where the dielectric constant changes) and the location of the jump ∆µ ± in solvation potential felt by the ions. When x < 0 (resp. x > 0) corresponds to the oil (resp. water) phase, assuming s > 0 and setting ψ(∞) = 0, we can write their result for the electrostatic potential as :

ψ P B =                        ψ D - ψ D A e κox , x < 0 ψ D - ψ D A [cosh(κ i x) + n sinh(κ i x)] ,
x

∈ [0, s] ψ D A e -κw(x-s) p[n cosh(κ i s) + sinh(κ i s)] ,
x > s

where

A = (1 + np) cosh(κ i s) + (n + p) sinh(κ i s), n = o / w , p = ρ o /ρ w = exp(-β∆µ av /2
) and the screening length in the intermediate region is

κ -1 i = (2βe 2 ρ o / w ) -1/2
. This analytical result is exact only in the linearized regime, i.e. potential differences small compared to k B T /e. Here we report lattice simulations based on our free energy model and compare the resulting ionic profiles to the prediction of the PB treatment for a flat interface.

System

We peformed simulations of flat oil-water interfaces in a simulation box consisting of N × 1 × 1 lattice points with N = 500 and containing two interfaces (one O/W and one W/O). Owing to the periodic boundary conditions, the system corresponds to an infinite stack of oil and water slabs of width 2L o and 2L w (we used L o = L w = 125∆x). There are several relevant lengthscales in the system: The interface width ξ, the Debye screening lengths κ -1 o,w and the size of both phases L o,w . The width of the interface should be small compared to all other lengthscales : κ o,w ξ 1 and ξ L o,w . In principle, the electric double-layers can overlap, especially in the oil phase where κ -1 o is larger. This is even a crucial point to explain the electrostatic stabilization of surfactant free water droplets in oil [START_REF] Leunissen | Proc. Natl. Acad. Sci. USA[END_REF][START_REF] De Graaf | [END_REF]. In the present paper, we want to assess the validity of our coarse-grained simulation method by comparison with known results for an isolated interface and therefore consider only the limit κ o,w L o,w → ∞. In particular, this condition is necessary for the salt concentrations to reach bulk values ρ o,w corresponding to well-defined screening lengths. Otherwise the amount of ions in the double layers could be non-negligible compared to the amount of ions in the bulk oil and water phases and the latter couldn't be considered as reservoirs. All simulations were performed using mixing free energy parameters βB∆x 3 = 10 -2 and βK∆x = 3 10 -2 with a mobility such that M B = 5 10 -2 ∆x 2 /∆t. These values of B and K lead in the absence of ions to an interface width ξ ∼ 2.45∆x. Results are given only for one of the two interfaces in the simulation box, with x < 0 (resp. x > 0) corresponding to the oil (resp. water) phase.

Results and discussion

FIG. 6: Ionic profiles at an oil (left) -water (right) interface, for γ = 0, β∆µ± = ±4, lB = 0.8∆x and ρ w ∆x 3 = 4.7 10 -4 . Symbols are simulation results. Lines are analytical results based on Eq. ( 27) with a shift s = 0 for the electrostatic potential ψ P B . Dashed lines corresponds to the fully lineared solution, solid lines to the reexponentiated one (see text). The insert compares the simulation results (symbols) for φ and ψ to the analytical solutions (lines).

We first investigate the influence of the solvation free energy differences ∆µ ± on the structure of the interface by considering a system without image charges ( w = o , i.e. γ = 0). Fig. 6 compares simulation results to analytical predictions based on Eq. ( 27) for β∆µ ± = ±4, corresponding to hydrophilic cations and hydrophobic anions. From the PB ionic chemical potential µ P B ± = ±eψ P B + V solv ± we compute the ionic concentrations either as ρ ± = ρ w e -βµ P B ± or as the linear expansion ρ w 1 -βµ P B ± . In the following, we refer to the former as the reexponentiated PB result (RPB, see e.g. Ref. [START_REF] Hansen | Theory of Simple Liquids[END_REF]) and to the latter as the linearized PB result (LPB). Since β∆µ av = 0 we expect equal salt concentrations in both phases and equal screening lengths (since we also have γ = 0) κ -1 o,w ∼ 10.3∆x, with a relatively good separation of length scales (κ o,w ξ < 1 and κ o,w L o,w > 1). We also expect a Donnan potential difference ψ D = -4k B T /e across the interface.

This figure shows that the LPB approximation fails to reproduce the ionic profiles, as expected for this large value of βeψ D . It even predicts negative concentrations near the interface. Adding the quadratic term in βeψ D to approximate ρ ± significantly improves the agreement while remaining consistent with the linear approximation for the potential (not shown). The RPB result, although not fully consistent with this approximation, is in quantitative agreement with the simulations, except at the interface where the finite width ξ smoothens the jump in ionic concentrations. The insert of Fig. 6 shows that the composition profile is close to φ th (x) = -tanh(x/ξ) and that the electrostatic potential is very well described by the linearized PB solution ψ P B (x) with a shift s = 0. Note that Eq. ( 27) for ψ P B is obtained under the assumption ξ = 0, and that potential differences are small compared to k B T /e. Surprizingly, simulations for a finite ξ with relatively large values of βeψ D give similar results for the electrostatic potential even in the transition region where the ionic profiles differ. Although it is not easily seen on this insert, the simulated profile is in fact sharper than φ th (x). Increasing βeψ D while keeping β∆µ av = 0 increases the charge separation and the excess of ions at the interface compared to the bulk concentration. Before analyzing quantitatively the influence of ψ D on these quantities, let us examine the effect of the finite interface width on the ionic profiles ρ ± (x) and the influence of the solvation free energy on the composition profile φ(x). Fig. 7 reports the ionic profiles for γ = 0.0, β∆µ av = 0 and Donnan potentials βeψ D = -1, -2 and -3. Simulation results are compared to the RPB solution for a thin interface (ξ = 0) and for solvation potentials V solv ± corresponding to a composition φ th (x) = -tanh(x/ξ) expected for an ion-free interface, where ξ = 2K/B. The main effect of the finite interface width is to smoothen ionic concentration profiles over a distance ∼ ξ. Beyond this distance the concentrations are equal to the thin interface result. The profiles are better described by the combination of ψ P B with the smoothed solvation potentials, but the agreement deteriorates as |ψ D | increases. This is because the composition profile φ(x) also deviates from φ th (x) as |ψ D | increases : The larger the solvation free energy difference (i.e. the larger |ψ D |), the sharper the interface.

We now analyze how the finite width ξ of the interface affects the overall excess of ions near the interface, quantified by the adsorption :

Γ = 0 -∞ [ρ + (x) + ρ -(x) -2ρ o ] dx + ∞ 0 [ρ + (x) + ρ -(x) -2ρ w ] dx . (28)
Fig. 8 compares the simulated adsorption to the LPB result :

Γ LP B = κ 8πl B (βeψ D ) 2 4 ( 29
)
with κ = κ o = κ w and l B the Bjerrum length common to both phases (since γ = 0). This result is obtained by expanding ρ ± (x) to second order in βµ P B ± , for the linear expansion yields Γ = 0. The error bars reported in Fig. 8 correspond to the estimate :

δΓ = 1 2 0 -∆x [ρ + (x) + ρ -(x) -2ρ o ] dx + ∆x 0 [ρ + (x) + ρ -(x) -2ρ w ] dx (30)
with ∆x the lattice spacing. The agreement between the simulation and LPB results is seen to be very good, except at very large ψ D . This performance of the LPB approximation even for βeψ D > 1 where it fails to predict the ionic profiles can be traced back to a compensation of errors. The RPB result, much closer to the simulated ionic profiles, predicts larger adsorptions than the LPB approximation. But the finite interface width smoothens the ionic profiles and therefore diminishes the value of Γ. In addition, for the largest ψ D , Γ is not small compared to the amount of "bulk" ions ρ o L o + ρ w L w ∼ 0.12∆x -2 so that the simulated system cannot be considered as an isolated interface in contact with infinite reservoirs. The good agreement between the simulated ψ and linearized PB result ψ P B suggests that the electric charge density is also close to the PB solution. In particular, the total electric charge of each phase σ el = σ o = -σ w = 0 -∞ (ρ + -ρ -)e dx can be compared to the PB prediction for ξ = 0 :

σ LP B el e = κ 4πl B A (βeψ D ) (31)
with here A = 2 (see Eq. ( 27)). 31) and confirms that the charge separation at the interface is proportional to ψ D , except at very large ψ D . The symmetric interface behaves as a capacitor of permittivity and width κ -1 , with capacitance per unit area κ. When submitted to a potential difference ψ D , each side builds up a charge per unit area σ el ∝ κψ D . The good agreement between the simulation and PB results for σ el is consistent with the finding of ref. [35] that global quantities related to the partitioning between the two phases are not in-fluenced by the finite interface width if it is smaller than the interfacial Debye length.

All results presented so far concerned the symmetric solvation case β∆µ av = 0. Simulations in the asymmetric case also give the expected results. For example, we find that the salt concentration ratio ρ o /ρ w decreases from e -1 ∼ 0.37 for β∆µ + = +2 and β∆µ -= 0 (hence β∆µ av = +1) to e -3 ∼ 5 10 -2 for β∆µ + = +4 and β∆µ -= +2 (β∆µ av = +3). While these two conditions correspond to the same ψ D = -k B T /e, the former yields a larger adsorption Γ. In addition, the ionic profiles are again in quantitative agreement with the RPB solution for a shift s = 0. This choice gives the best agreement when there is no dielectric contrast between the two solvents.

FIG. 9: Ionic profiles at the oil (left) -water (right) interface for γ = 0.5, β∆µ± = ±2, lB = 0.8∆x, ρ w ∆x 3 = 4.9 10 -4 and a shift s = -1.25∆x. Simulation results (symbols) are compared to analytical results for a thin interface (solid line) and for the combination of ψ P B with solvation potentials V solv ± corresponding to a composition φ th (x) = -tanh(x/ξ) (dashed line). The insert compares the simulation (symbols) and analytical (lines) results for the composition φ and electrostatic potential ψ.

When γ = 0 the situation is more complex and even the sign of s depends on γ, ∆µ + and ∆µ -. An example of such a situation is illustrated in Fig. 9. The best agreement with Eq. ( 27) was obtained in that case for a shift s = -1.25∆x. All other parameters being fixed, we find that for increasing ∆µ + = -∆µ - the required shift also increases. Simulation results for the electrostatic potential are well described by the PB result for a vanishing width ξ = 0 and a finite shift. Therefore the minimal model of van Roij et al.. to account for image charges seems to be appropriate, at least for the explored range of parameters. Even for βeψ D > 1, where the fully linear PB approximation fails, we find that the reexponentiated PB result gives an accurate description of ionic profiles. Moreover, our results show that the finite width of the interface affects the ionic profiles in such a way that the electrostatic potential is still close to the PB prediction. These profiles are accurately described by the combination of ψ P B with solvation potentials V solv ± corresponding to an unperturbed composition profile φ th (x) = -tanh(x/ξ). In particular, "far" from the interface, the effect of the finite width on solvation forces is negligible and the ionic profiles coincide with the result for ξ = 0. As shown previously in the case γ = 0, we expect for larger ψ D a perturbation of the composition profile and a corresponding modification of the ionic concentrations.

The results presented in this section demonstrate the ability of our coarse-grained simulation scheme to study the interface between two immiscible solvents in presence of ions, including the effect of their (possibly asymmetric) affinity for both phases. To assess the validity of the simulation scheme we only showed results in the one-dimensional case, for which analytical predictions are available. Simulation in two or three dimensions is straightforward.

CONCLUSION

We have presented a coarse-grained simulation method for complex charged systems. This mesoscopic model couples a hydrodynamic description to a free energy functional accounting for the interactions between solvent(s) and charged solutes. All the parameters entering in the model, such as free energy parameters (e.g. related to ion solvation), solvent viscosity or ionic diffusion coefficients are, at least in principle, computable by simulations at the molecular level. We described the implementation of this model in a hybrid lattice-based scheme, whereby the evolution of the overall mass and momentum is taken care of via a Lattice Boltzmann scheme, whereas the composition and ionic concentrations are updated using the linkflux method. We presented several applications of the coarse-grained simulation method: the transport of charged tracers in charged porous media, the deformation of an oil droplet in water under the effect of an applied electric field, and the distribution of ions at an oil-water interface as a function of the ions affinity for both solvents. When possible, we compared our simulation results to exact or approximate analytical results to investigate the range of simulation parameters that can be used to recover the continuous results (e.g. the effect of a finite interface width).

The proposed method will be very useful to sim-ulate the dynamics of complex mixtures of solvents and ions. In particular, it will be interesting to investigate electrokinetic phenomena at a charged oilwater interface, or of oil-water mixtures in charged porous media. Understanding electrokinetic effects in these systems might be very helpful in designing new electro-acoustic oil recovery techniques or monitoring devices. This coarse-grained simulation method could also be applied to the study of electrowetting and of microfluidic devices, particularly those based on electrokinetic pumping.

Although surfactant free emulsions can be stabilized by the presence of salt, the most usual situation also involves sufactant molecules at the interface and in solution. As free energy based models for mixtures of oil, water and surfactants have already been proposed, it should be rather straightforward to couple them to the one introduced in the present paper. One physical feature not included in the method presented here is the presence of thermal fluctuations. It would be interesting to investigate the possibility of including fluctuations as is done in the fluctuating lattice-Boltzmann method [START_REF] Adhikari | [END_REF]64].

2πlB 1 cos 2

 12 αx , where α satisfies αL 2 tan αL 2 = πσ el Ll B /e. If an electric field of strength E 0 is applied parallel to the solid walls, an electroosmotic flow develops in the confined fluid, u y (x) = u ref ln cos αx cos αL/2 (22)

FIG. 1 :

 1 FIG.1: Steady state flow profile across a slit of width L = 60∆x (with lB = 0.4∆x) for an applied electric field of magnitude βeE0L = 3 as a function of the the solid surface charge density. Symbols are simulation results while the continuous curves correspond to the theoretical prediction.

2 ∼

 2 tion (porosity) of 1 -π/3 √ 26% is divided into large octahedral (O h ) cavities of radius r O h ∼ 0.41R connected by smaller tetrahedral (T d ) pores of radius r T d ∼ 0.22R. The size of the bottlenecks between O h and T d pores is approximately 0.15R. The electrostatic potential distribution inside the pore is controlled by the salt concentration, with a typical double-layer thickness κ -1 b = (8πl B ρ b ) -1/2 , with l B the Bjerrum length. The Debye length κ -1 b

FIG. 2 :

 2 FIG. 2: Effective diffusion coefficient for charged tracers, as a function of the salt concentration. The results for neutral (•), cationic ( ) and anionic ( ) tracers , normalized by the molecular diffusion coefficient D, are reported as a function of the equivalent Debye length in a bulk solution κ -1 b = (8πlBρ b ) -1/2 divided by the radius R of the spheres. The insert shows the characteristic time to explore the porosity accessible to each tracer, normalized by R 2 /D.

  with σ the surface tension and R d the drop radius. When the drop and the suspending liquid have different dielectric constants (i.e. when γ = wo w + o = 0), applying an electric field E polarizes the drop and the anisotropic electrostatic stress tensor tends to elongate it in the direction of the field.

FIG. 4 :

 4 FIG. 4: Electrostatic potential (background) and electric field lines (the color of streamlines indicates the magnitude of the field, with red for large fields) for an oil droplet in water under an electric field. The oil-water interface is indicated as a dashed line. The dielectric constrast is γ = 0.9. Analytical results for an isolet droplet (a and c) are compared to simulation results with periodic boundary conditions (b and d). The effect of the finite width (ξ = 2.45∆x for b and d) is more visible for the smaller system (a,b: box size L = 50∆x, droplet radius R d = 14.1∆x) than with the larger one (c,d: L = 100∆x, R d = 28.2∆x).

FIG. 5 :

 5 FIG. 5: Deformation of a 2D-oil droplet in water, as a function of the dielectric contrast γ = ( wo)/( w + o). The deformation D is normalized by ¯ E 2 R d /σ, with ¯ = ( w + o)/2, E the applied field, R d the droplet radius and σ the surface tension. Simulation results are for two droplet radii (R d = 14.1∆x for , × and 28.2∆x for •, +) and two electric field strengths (βeE∆x = 10 -3 for •, and 3 10 -3 for ×, +) with the same interface width ξ = 2.45∆x. The solid line corresponds to the analytical solution (25).

FIG. 7 :

 7 FIG. 7: Ionic profiles in the water phase for γ = 0.0, β∆µav = 0, lB = 0.8∆x, ρ w ∆x 3 = 4.9 10 -4 and Donnan potentials βeψD = -1 (a), -2 (b) and -3 (c). Simulation results (symbols) are compared to the RPB result (see text) for a thin interface (solid line) and for V solv ± corresponding to φ th (x) = -tanh(x/ξ) (dashed line) expected for an ion-free interface.

FIG. 8 :

 8 FIG. 8: Ionic adsorption Γ and surface charge density σ el of the oil phase as a function of ψD for γ = 0 and β∆µav = 0. Simulation results (symbols) for lB = 0.8∆x are compared to the linearized Poisson-Boltzmann results for ξ = 0 (lines) given by Eqs. (29) and (31).

Fig. 8 2 ∆x-

 82 also displays σ el as a function of the Donnan potential ψ D for the simple case β∆µ av = 0 and γ = 0. The error bars correspond to δσ el = 1 ∆x |ρ + -ρ -| e dx. This figure shows that simulation results are indeed well described by Eq. (

Acknowledgements

The authors would like to thank Rene van Roij, Jean-Pierre Hansen and Edo Boek for very useful discussions. B.R. acknowledges financial support from the Agence Nationale pour la Gestion des Déchets Radioactifs (ANDRA, France). I.P. acknowledges financial support from CAICYT (FIS2005-01299). The work of the FOM Institute is part of the research program of FOM and is made possible by financial support from the Netherlands organization for Scientific Research (NWO).